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Abstract 

Benthic  shells  can  contribute  greatly  to  the  scattering  variability  of  the  ocean  bottom, 
particularly  at  low  grazing  angles.  Among  the  effects  of  shell  aggregates  are  increased 
scattering  strength  and  potential  subcritical  angle  penetration  of  the  seafloor.  Sand  dollars 
(. Dendraster  excentricus )  occur  commonly  in  the  ocean  and  have  been  shown  to  be  significant 
scatters  of  sound.  In  order  to  understand  more  fully  the  scattering  mechanisms  of  these 
organisms,  the  scattering  from  individual  sand  dollars  was  studied  using  several  methods. 

Using  an  approximation  to  the  Helmholtz-Kirchhoff  integral,  the  Kirchhoff  method 
gives  an  analytic  integral  expression  to  the  backscattering  from  an  object.  This  integral  was  first 
solved  analytically  for  a  disk  and  a  spherical  cap,  two  high  aspect  ratio  oblate  shapes  which 
simplify  the  shape  of  an  individual  sand  dollar.  A  method  for  solving  the  Kirchhoff  integral 
numerically  was  then  developed.  An  exact  three  dimensional  model  of  a  sand  dollar  test  was 
created  from  computed  tomography  scans.  The  Kirchhoff  integral  was  then  solved  numerically 
for  this  model  of  the  sand  dollar. 

The  finite  element  method,  a  numerical  technique  for  approximating  the  solutions  to 
partial  differential  equations  and  integral  equations,  was  used  to  model  the  scattering  from  an 
individual  sand  dollar  as  well.  COMSOL  Multiphysics  was  used  for  the  implementation  of  the 
finite  element  method. 

Modeling  results  were  compared  with  published  laboratory  experimental  data  from  the 
free  field  scattering  of  both  an  aluminum  disk  and  a  sand  dollar.  Insight  on  the  scattering 
mechanisms  of  individual  sand  dollar,  including  elastic  behavior  and  diffraction  effects,  was 
gained  from  these  comparisons. 
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Chapter  1  Introduction 


1.1  Sound  as  a  Tool  in  the  Ocean 

Sound  is  an  important  tool  for  gaining  insight  into  phenomena  and  processes  in 
underwater  environments.  In  the  ocean,  the  theory  of  sound  is  essential  in  such  broad 
applications  as  military  operations,  geological  exploration,  and  biological  surveys.  Unlike 
electromagnetic  waves,  acoustic  pressure  waves  are  able  to  propagate  long  distances  in  water. 
The  typical  frequency  range  used  in  underwater  acoustics  is  10  Hz  to  1  MHz,  with  the  lower 
frequencies  able  to  travel  many  kilometers.  An  important  aspect  of  the  field  is  acoustic 
scattering,  that  is,  understanding  how  sound  reacts  to  boundaries  and  obstacles  in  its  path.  With 
SONAR  (Sound  navigation  and  ranging)  —  a  method  for  marine  vessels  to  navigate, 
communicate,  and  detect  one  another  — as  its  most  well-known  application,  acoustic  scattering 
includes  the  study  of  the  reflection,  diffraction,  and  transmission  of  sound  incident  upon  a 
particular  object.  This  analysis  can  often  convey  detailed  information  about  the  nature  of  the 
object  such  as  its  shape,  size,  or  material  properties.  Knowledge  of  scattering  mechanisms  is 
important  in  such  diverse  applications  as  mine  detection  and  investigating  zooplankton 
populations.  There  is  a  vast  literature  on  the  subject  of  underwater  acoustic  scattering  (Urick, 
1983;  Medwin  and  Clay,  1998). 
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1.2  Scattering  from  Benthic  Shells 


The  seafloor  is  a  significant  scatterer  of  sound.  Not  only  essential  to  the  determination  of 
bathymetry,  knowledge  of  how  sound  scatterers  from  the  ocean  bottom  has  commercial, 
military,  and  scientific  ramifications.  Sediment  type,  bottom  roughness,  and  grazing  angle  all 
are  important  variables  in  seafloor  scattering.  Jackson  et  al.  (1986)  presented  data  showing  how 
sediment  type  alone  is  not  a  good  predictor  of  scattering  strength,  with  variations  up  to  15  dB 
within  a  single  class  of  seafloor  material.  Jackson  and  colleagues  also  showed  how  seafloor 
roughness  seemed  to  account  for  high  variability  in  the  acoustic  returns  at  shallow  grazing 
angles.  An  important  instance  of  seafloor  roughness  is  the  presence  of  benthic  shelled 
organisms.  Studies  have  suggested  that  the  presence  of  benthic  shells  can  contribute  greatly  to 
the  scattering  strength  of  the  ocean  bottom  at  frequencies  above  10  kHz  (Stanic  et  al,  1989; 
Fenstermacher  et  al,  2001).  Not  only  can  shells  greatly  increase  the  scattering  strength  in  the 
backwards  direction,  but  they  also  scatter  sound  into  the  seafloor  itself.  This  raises  the 
possibility  of  detecting  buried  objects  at  subcritical  angles  (Stanton  and  Chu,  2004).  However, 
the  understanding  of  benthic  shell  aggregate  scattering  remains  limited  due  to  an  incomplete 
knowledge  of  the  scattering  from  individual  benthic  shells. 

Accurately  describing  the  scattering  from  most  marine  organisms  is  a  significant 
challenge,  in  part  due  to  their  geometric  complexities.  In  particular,  there  are  no  exact  analytic 
expressions  for  the  scattering  from  irregular  benthic  shells.  There  are  only  eleven  relatively 
simple  geometric  shapes  with  separable  coordinate  systems,  including  spheres  and  infinite 
cylinders,  and  thus  with  exact  analytic  solutions  to  their  scattering  (Bowman  et  al,  1969; 
Anderson,  1950;  Faran,  1951).  In  addition  to  geometry,  material  properties  are  also  an 
important  factor  in  the  scattering  physics  of  an  object.  These  material  properties  are  important 
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in  determining  proper  boundary  conditions  on  the  surface  of  the  object.  Fluid  objects  simply 
require  continuity  of  pressure  and  velocity  on  their  boundary.  Solid  objects,  referred  to  as 
elastic,  involve  more  complicated  boundary  conditions  due  to  the  presence  of  shear  waves. 
Accurately  representing  the  boundary  conditions  of  an  object  can  be  just  as  important  as  its 
geometry  in  the  modeling  of  its  scattering. 

There  has  been  much  work  done  regarding  the  scattering  of  individual  zooplankton, 
including  shelled  organisms,  in  the  water  column  (Simmonds  and  MacLennan,  2005).  Certain 
animal  classes  have  been  approximated  as  simple  geometrical  shapes,  such  as  spheres 
(Greenlaw,  1977;  Johnson,  1977;  Stanton  et  al,  1987;  2000)  and  finite  cylinders  (Stanton,  1989; 
Stanton  et  al.r  1993).  More  complex  representations  of  certain  water-column  scatterers  have  also 
been  used.  In  particular,  scattering  models  based  on  the  distorted  wave  Born  approximation 
(DWBA)  have  proved  useful  in  modeling  fluid-like  zooplankton  (Chu  et  al.,  1993;  Stanton  et  al., 
1993, 1998;  Stanton  and  Chu,  2000;  Lavery  et  al,  2002;  Lawson  et  al,  2006,  Lavery  et  al,  2007). 

Compared  to  the  understanding  developed  for  scattering  from  water-column  organisms, 
relatively  little  work  has  been  done  regarding  the  scattering  from  individual  benthic  shelled 
organisms.  Stanton  et  al.  (2000)  used  a  deformed  sphere  in  both  ray-based  and  modal-series- 
based  approaches  to  model  scattering  from  periwinkles  ( Littorina  littorea),  a  type  of  benthic 
shelled  animal.  Stanton  and  Chu  (2004)  compared  laboratory  measurements  of  scattering  from  a 
machined  round,  elastic  disk  to  that  of  a  sand  dollar  ( Dendraster  excentricus)  and  bivalve 
(. Dinocardium  robustum  vanhyningi).  Specifically,  they  observed  diffraction  effects  from  the  edges 
of  the  sand  dollar  through  pulse-compression  analysis.  Though  much  progress  was  made  in 
these  studies,  understanding  the  importance  of  the  geometrical  complexities  of  these  benthic 
shelled  organisms  on  the  scattering  was  incomplete,  which  is  the  motivation  for  this  study.  In 
particular,  the  focus  of  this  work  is  to  understand  acoustic  scattering  from  sand  dollars,  for 


13 


which  there  are  controlled  laboratory  scattering  data  available  (Stanton  and  Chu,  2004)  that  can 
give  insight  and  inspire  model  development. 


1.3  Scattering  from  Sand  Dollars 

Sand  dollars  ( Dendraster  excentricus)  were  shown  by  Fenstermacher  et  al.  (2001)  to 
contribute  significantly  to  seafloor  scattering  levels  when  in  naturally  occurring  dense 
collections.  These  benthic  echinoderms  can  form  concentrations  of  up  to  several  hundred  per 
square  meter  in  the  sandy,  shallow  water  of  the  west  coast  of  North  America  (Chia,  1969; 
Highsmith,  1982).  Sand  dollars  have  skeletons  composed  of  strong  calcite,  correctly  referred  to 
as  tests  and  not  shells  because  of  a  small  amount  of  external  living  tissue  (Nichols,  1969).  These 
tests,  clearly  visible  after  the  animal  has  died,  resemble  thin  round  disks,  as  seen  in  figure  1. 


Figure  1  -  Photograph  of  sand  dollar  test,  top-view  (Stanton  and  Chu,  2004). 

Sand  dollars  grow  to  several  centimeters  in  diameter;  Chia's  (1969)  population  records  show  an 
average  adult  diameter  of  7  cm.  Understanding  the  scattering  from  individual  sand  dollars  is 
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important  for  predicting  the  scattering  from  aggregations  of  live  organisms  found  in  situ  on  the 
ocean  bottom. 

This  work  attempts  to  further  understand  the  scattering  from  individual  sand  dollars  by 
first  modeling  them  as  high  aspect  ratio  oblate  objects,  namely  disks  and  spherical  caps.  Along 
with  the  geometrical  complexities,  the  material  properties  of  sand  dollars  are  important  to  an 
accurate  description  of  their  scattering.  One  key  simplification  to  the  models  developed  here  is 
that  a  rigid/ fixed,  rather  than  elastic,  boundary  condition  is  used.  This  removes  the 
complexities  of  shear  waves  and  internal  structure  from  the  analysis  and  proves  to  be  a  good 
approximation  at  certain  angles  of  orientation  with  an  added  heuristic  reflection  coefficient. 
Though  the  complex  scattering  mechanisms  of  elastic  disks  have  been  illustrated  by  Hefner  and 
Marston  (2001),  they  are  beyond  the  scope  of  this  work.  Also,  the  models'  focus  is  solely  on  the 
backscattering  of  incident  harmonic  plane  waves.  Both  an  approximate  analytical  technique, 
namely  the  Kirchhoff  method,  and  a  fully  numerical  technique,  namely  the  finite  element 
method,  are  used  to  model  the  scattering  from  sand  dollars.  Both  techniques  model  them  as 
simple  high  aspect  ratio  oblate  objects  (disks  and  spherical  caps).  The  Kirchhoff  method  is  then 
taken  further  to  model  the  scattering  from  the  sand  dollar's  exact  geometry,  obtained  from  high 
computed  tomography  (CT)  scans.  All  work  involves  free  field  scattering,  that  is,  the  scatterer  is 
away  from  all  boundaries. 


1.4  The  Kirchhoff  Method 

The  Kirchhoff  method  is  an  approximation  of  the  Helmholtz-Kirchhoff  integral,  derived 
from  the  theorems  of  Gauss  and  Green.  The  result  of  the  method  is  an  analytic  integral  that  is 
evaluated  over  the  insonified  surface  of  the  scatterer  to  determine  its  reflection.  There  is  a  vast 
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literature  on  the  subject  in  both  optics  and  acoustics  (Born  and  Wolf,  1970;  Medwin  and  Clay, 
1998).  Because  any  scattering  process  often  also  involves  strong  diffraction  effects,  there  are 
limitations  to  the  accuracy  of  the  method.  In  particular,  it  loses  accuracy  as  wavelength 
increases  beyond  the  characteristic  dimensions  of  the  object.  Norton  et  al.  (1993)  studied  the 
diffraction-induced  errors  of  this  method  in  the  scattering  from  circular  disks. 

In  Chapter  2,  a  brief  derivation  of  the  Kirchhoff  method  is  discussed.  The  resulting 
integral  expression  is  then  solved  for  two  objects:  a  sphere  and  a  finite  cylinder.  The  results  are 
compared  to  the  sphere's  exact  modal  series  solution  and  the  finite  cylinder's  modal-series- 
based  approximation.  Next,  the  Kirchhoff  integral  expression  is  solved  for  a  disk  and  spherical 
cap,  the  two  objects  that  approximate  the  shape  of  each  side  of  the  sand  dollar.  Finally,  a 
method  for  solving  the  Kirchhoff  method's  integral  numerically  is  presented.  It  is  based  on 
creating  a  mesh  conforming  to  the  surface  of  the  scatterer  and  summing  the  integral  over  the 
resulting  elements.  This  approach  is  valuable  because  it  does  not  require  the  object's  surface  to 
be  separable  and  thus  can  be  used  on  arbitrarily  complex  scatterers.  Finally,  the  numerical 
approach  is  tested  against  the  previous  analytic  solutions  of  the  Kirchhoff  integral  for  the 
sphere,  finite  cylinder,  disk,  and  spherical  cap. 


1.5  The  Finite  Element  Method 

Broadly  defined,  a  numerical  technique  is  a  way  of  approximating  continuous  problems 
with  discrete  solutions.  The  idea  of  numerical  analysis  arose  long  before  computers,  evidenced 
by  very  old  techniques  such  as  Newton's  and  Euler's  methods.  However,  the  surge  in 
computing  power  over  the  past  several  decades  has  greatly  expanded  the  scope  and  power  of 
using  numerical  techniques  for  solving  problems.  In  the  context  of  acoustic  scattering,  a  few 
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distinct  numerical  techniques  have  risen  in  popularity,  such  as  the  T-matrix  and  finite  element 
methods.  The  T-matrix  method  is  a  formally  exact  numerical  solution  to  the  wave  equation. 
Introduced  by  Waterman  (1969),  this  approach  has  been  used  to  study  the  scattering  and 
diffraction  of  thin  circular  disks  (Kristensson  and  Waterman,  1982;  Norton  et  al.,  1993;  Stanton  et 
al,  2007). 

An  approach  to  approximating  solutions  to  partial  differential  equations  and  integral 
equations,  the  finite  element  method  is  a  general  numerical  technique  with  a  broad  scope  of 
applications.  It  stemmed  from  the  work  of  applied  mathematicians,  physicians,  and  engineers 
in  the  mid  twentieth  century;  each  community  was  interested  in  approximating  solutions  to 
increasingly  complex  continuous  problems  (Huebner,  1975).  The  first  mention  of  the  term 
"finite  element  method"  was  by  Clough  (1960)  to  describe  plane  elasticity.  Earlier,  Courant 
(1943)  used  piecewise  continuous  functions  over  triangular  elements  to  study  the  St.  Venant 
torsion  problem.  The  mathematical  foundations  of  the  finite  element  method  were  solidified  in 
the  early  1970s  (Strang  and  Fix,  1973;  Babuska  and  Aziz,  1973),  and  it  has  become  a  popular 
numerical  technique  in  a  wide  variety  of  fields.  Early  use  of  the  finite  element  method  in  the 
study  of  acoustics  focused  on  internal  problems  such  as  acoustics  of  an  enclosure,  structural 
vibration,  and  waveguides  (Gladwell,  1966;  Craggs,  1972;  Nefske  et  al,  1982;  Petyt,  et  al,  1976). 
As  techniques  have  been  developed  to  deal  with  the  infinite  domain  of  scattering  problems,  the 
finite  element  method  has  been  adapted  for  them  as  well  (Hunt  et  al.,  1975;  Bettess,  1977;  Gan  et 
al,  1993;  Berenger,  1994;  Thompson,  2006). 

In  Chapter  3,  the  theory  behind  the  finite  element  method  and  its  application  to  acoustic 
scattering  problems  are  discussed.  Techniques  for  reducing  numerical  error  in  two  and  three 
dimensions  are  presented.  The  program  COMSOL  Multiphysics®  is  then  explored  for  its  use  in 
predicting  the  scattering  from  sand  dollars.  The  accuracy  and  limitations  of  this  method  are 
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explored  by  again  solving  problems  with  known  solutions,  scattering  from  a  rigid/fixed  infinite 
cylinder  and  a  rigid/fixed  sphere. 


1.6  Comparison  to  Experiment 

In  Chapter  4,  the  results  from  both  the  Kirchhoff  method  and  finite  element  method 
models  are  compared  with  published  (Stanton  and  Chu,  2004)  experimental  results.  In  Stanton 
and  Chu  (2004),  forward  scattering  and  backscattering  from  a  sand  dollar  test,  a  bivalve  shell, 
and  a  machined  aluminum  disk  of  similar  size  were  measured  over  a  range  of  angles  of 
orientation.  Model  results  based  on  the  finite  element  method  and  analytic  and  numerical 
solutions  to  the  Kirchhoff  method  are  compared  to  the  experimental  results  from  the  aluminum 
disk  and  the  two  faces  of  the  sand  dollar.  These  comparisons  are  conducted  at  a  single 
frequency,  70  kHz,  over  a  range  of  angles  of  incidence. 

As  the  wavelength  used  is  similar  in  size  to  the  dimensions  of  the  sand  dollar  and 
aluminum  disk,  certain  inaccuracies  show  up  in  the  use  of  the  Kirchhoff  method,  as  mentioned 
before.  However,  near  normal  incidence  to  the  main  faces  of  the  sand  dollar  and  aluminum  disk 
(referred  to  as  broadside),  this  approach  results  in  good  agreement  between  the  predicted  and 
measured  scattering.  The  finite  element  method  is  able  to  model  some  of  the  experimental  data 
at  higher  angles  of  orientation  in  addition  to  working  well  near  broadside.  Errors  due  to  the  use 
of  a  rigid/fixed  boundary  condition  are  discussed  and  a  heuristic  approach  to  representing  the 
true  penetrable  condition  of  the  scatterers  by  means  of  an  added  reflection  coefficient  is 
developed. 
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1.7  Overview  of  Study 


To  summarize.  Chapter  2  discusses  the  Kirchhoff  method  and  the  resulting  solutions  for 
two  shapes  similar  to  a  sand  dollar.  Chapter  3  discusses  the  finite  element  method  and  its 
application  in  modeling  the  scattering  by  a  sand  dollar.  Chapter  4  compares  the  results  of  the 
models  with  experimental  scattering  data  for  an  aluminum  disk  and  sand  dollar.  A  heuristic 
correction  to  the  models  is  discussed  as  well  to  account  for  the  use  of  the  rigid / fixed  rather  than 
elastic  boundary  conditions.  Chapter  5  summarizes  the  thesis  and  offers  recommendations  for 
future  work. 
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Chapter  2  The  Kirchhoff  Method 


2.1  Introduction 

In  this  chapter,  models  for  understanding  the  backscattering  from  individual  sand  dollars 
are  created  through  use  of  the  Kirchhoff  method.  A  brief  overview  of  the  Kirchhoff  method  is 
presented,  which  results  in  an  analytic  integral,  referred  to  as  the  Kirchhoff  integral,  which 
approximates  the  backscattering  amplitude  of  an  object  when  exposed  to  an  incident  harmonic 
plane  wave.  The  Kirchhoff  integral  is  tested  for  two  shapes  with  known  solutions,  namely 
spheres  and  finite  cylinders.  However,  because  sand  dollars  have  complicated  shapes,  for 
which  there  are  no  analytic  solutions  to  the  Kirchhoff  integral,  they  are  represented  in  this 
chapter  by  objects  with  simpler  geometries  that  incorporate  the  sand  dollar's  high  aspect  ratio 
oblate  shape,  namely  round  disks  and  spherical  caps.  These  objects  are  all  rigid/fixed,  based  on 
the  assumption  that  this  is  a  good  approximation  of  the  sand  dollar's  boundary  condition. 
Finally,  a  technique  for  solving  the  integral  numerically  for  any  shape,  regardless  of  geometric 
complexities,  is  presented.  Numerical  solutions  to  the  Kirchhoff  integral  are  then  compared  to 
the  analytic  solutions  for  a  sphere,  finite  cylinder,  disk,  and  spherical  cap. 


2.2  Theoretical  Background 

The  Kirchhoff  method's  derivation  begins  with  the  theory  of  wave  propagation 
proposed  by  Christian  Huygens.  Described  in  detail  by  Born  and  Wolf  (1970),  Huygens' 
principle  is  that  each  point  on  an  advancing  wave  front  can  be  considered  a  source  of  secondary 


20 


spherical  wavelets.  Together,  these  wavelets  form  the  new  wave  front  as  it  propagates.  Using 
the  theorems  of  Gauss  and  Green,  analytic  expressions  for  these  wavelets  can  be  written.  In 
scattering  problems,  these  wavelets  are  analyzed  on  the  surface  of  the  scatterer,  resulting  in  the 
scattered  pressure  field.  Medwin  and  Clay  (1998)  summarize  the  derivation,  giving  the 
Helmholtz-Kirchhoff  integral. 


Pscat(x) 
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(2.1) 


This  is  an  expression  for  the  scattered  pressure  at  position  x  in  terms  of  its  value  along  the 
surface  of  the  scatterer,  o .  Each  differential  element  do  on  the  object's  surface  has  position  r' 
and  outward  unit  normal  h  .  The  term  ra  is  the  distance  from  the  surface  element  to  position  x , 
ra  =1 3c  — r  I  and  k  (—2ft  /  A ,  where  A  is  the  wavelength)  is  the  acoustic  wavenumber.  The 
Kirchhoff  approximation  simplifies  the  integral  by  assuming  that  each  differential  element 
reflects  sound  in  a  ray-like  manner;  each  element  behaves  just  as  a  part  of  a  tangent  infinite 
planar  interface  would  (Medwin  and  Clay,  1998).  This  means  that  the  scattered  pressure  can  be 
written  in  terms  of  the  incident  pressure  and  the  reflection  coefficient  Rj  that  would  result  from 
such  an  interface.  Implicit  in  this  approximation  is  that  only  the  insonified  or  "front"  of  the 
object  contributes  to  scattering.  The  field  in  the  shadow  or  "back"  of  the  object  is  zero. 

Therefore,  the  integral  for  the  scattered  pressure  is  only  evaluated  on  the  scatter's  insonified 
surface. 
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d/7 


(  JK  \ 

Pinc(r)  — 
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V  J 


To  simplify  further,  the  incident  field  is  considered  to  be  a  harmonic  plane  wave,  written  in  its 
time  independent  form  as  pinc(f)  =  p0e,kr ,  with  amplitude  pQ  and  wave  vector  k  ,  wire  re 
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k  =  kk  and  k  is  a  unit  vector  in  the  direction  of  the  wave's  propagation.  (Throughout  this 
work,  only  harmonic  waves  are  considered,  so  the  time  dependence  is  omitted).  Inserting  the 
expression  for  the  incident  field  into  the  integral  gives 


"'v  f  ikr ■  \ 

(3c)  =  -^-f  Rf  —  e*'  —  da. 
4n{  c)n 


Two  further  simplifications  are  made.  First,  the  scatterer  is  considered  to  have  a  rigid/ fixed 
boundary  condition;  this  is  defined  as  one  where  both  the  gradient  of  pressure  and  the  fluid 
velocity  are  zero  on  the  boundary.  The  reflection  coefficient  for  a  rigid/fixed  interface  is  always 
equal  to  one,  regardless  of  angle  of  orientation.  Second,  the  scattered  pressure  is  only  calculated 
at  a  distance  from  the  object  much  larger  than  the  object's  dimensions.  To  do  this,  the 
coordinate  system  must  be  oriented  so  that  the  scatterer  is  located  at  the  origin.  Then  the 
distance  from  the  origin  to  the  point  x ,  defined  as  r  =1  x  I,  must  satisfy  r»d,  where  d  is  a 
typical  object  dimension.  When  this  is  the  case,  r  is  a  good  approximation  of  the  ra  term  in  the 
denominator.  The  ra  in  the  exponent  must  be  approximated  more  carefully  due  to  the 
importance  of  the  phase, 

-  x  (2-4) 

ra~  r  —  r  ■  — . 

r 

When  the  scattered  pressure  is  calculated  in  the  backwards  direction,  opposite  to  the  direction 
of  the  incident  wave,  this  can  be  rewritten  as, 

ra  ~  r+r' -k  .  (2-5) 

With  these  approximations,  the  backscattered  pressure  (from  equation  (2.3))  can  be  expressed  as 
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The  normal  derivative  of - is  assumed  to  be  zero  because  r  is  large,  giving 
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The  integral  thus  reduces  to 


ikr 


Pbs(x)  =  jPo^y  J  (* '  k)e'2k"  da . 


(2.8) 


This  is  the  final  Kirchhoff  method  expression  for  the  backscattered  pressure. 
In  the  far-field  kr  »  1 ,  the  scattered  pressure  is  defined  as 


ikr 

Psca,  =  — 

r 


(2.9) 


where  f  (ki)  is  known  as  the  scattering  amplitude.  It  is  a  function  of  the  spherical  angles 
represented  by  T> .  Thus  from  equation  (2.8),  the  backscattering  amplitude,  fhs ,  can  be  written. 


fhs=j\(k-n)e'2krd(T. 


(2.10) 


Backscattering  amplitude  can  be  used  to  find  target  strength  (Urick,  1983), 

r^ioiogl/J2.  (zn> 

The  units  of  target  strength  are  decibels  (dB)  relative  to  1  m2.  Target  strength  is  a  heavily  used 
indicator  of  scattering  strength  and  is  based  solely  on  the  properties  of  the  scatterer  itself.  Only 
backscattering  is  considered  throughout  this  work,  so  the  subscript  bs  can  be  dropped.  This 
expression  in  equation  (2.10)  is  referred  to  here  as  the  Kirchhoff  integral  and  can  be  used  to 
calculate  the  scattering  of  an  object  by  performing  the  integral  over  its  insonified  surface. 

The  Kirchhoff  method  is  also  referred  to  as  the  physical  optics  method.  By  only 
accounting  for  reflection  from  the  object's  front  insonified  surface,  this  method  neglects  the 
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effects  of  diffraction.  Diffraction  effects  vary  greatly  with  frequency  and  object  shape,  making 
the  Kirchhoff  method  a  poor  approximation  at  times.  These  effects  can  be  explored  by 
comparisons  with  shapes  that  have  known  exact  solutions  or  approximations  to  their  scattering, 
such  as  spheres  and  finite  cylinders. 


2.3  Scattering  from  a  Sphere 

In  this  section,  the  Kirchhoff  integral  is  solved  for  a  rigid / fixed  sphere.  The  results  are 
compared  to  the  exact  modal  series  solution  for  scattering  from  a  rigid/fixed  sphere. 


2.3.1  Kirchhoff  Integral  Solution  for  a  Sphere 

The  first  step  involves  setting  up  a  spherical  coordinate  system  with  radial  variable  p , 
polar  variable  0 ,  and  azimuth  variable  (p .  The  Cartesian  coordinate  system  is  described  with 
unit  vectors  x ,  y ,  and  z  .  Assume  a  rigid/ fixed  sphere  with  radius  a  is  centered  at  the  origin 
and  consider  an  incident  plane  wave  pinc  =  p0e~lkz  proceeding  in  the  negative  z  direction,  as 
seen  in  figure  2.  This  orientation  is  chosen  to  simplify  the  limits  of  integration  in  the  integral. 


24 


z 


Figure  2  -  Spherical  coordinate  system,  shown  with  incoming  wave  and  without  scatterer  geometry. 


It  is  seen  that 


k  =  -kz 

h  =  cos^sin#x  +  sin^sin#y  +  cos  Oz 
r'  =  a  cos  (p sin  Ox  +  a  sin  (p  sin  6y  +  a  cos  Oz 
do  =  cc  smOdOdcp . 

In  this  case,  the  Kirchhoff  integral  becomes 

.  2  K  Kl  2 

f=-\  J  -cos Oe~i2kacosea 2 sin OdOdcp . 


0  0 


The  integral  over  (p  is  easily  performed,  leaving 
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(2.13) 


(2.14) 


Making  the  substitution  of  x  =  cos  0 , 
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f  =  -ika2\-xe~i2kaxdx 


(2.15) 


the  expression  for  backscattering  amplitude  becomes 


•  -i2ka  •  -ilka 

r  i  ae  le 

i = — + - 

4k  2  4  k 


(2.16) 


Target  strength  follows  from  equation  (2.11). 


2.3.2  Modal  Series  Solution  for  a  Sphere 

Again  consider  the  spherical  coordinate  system  with  a  rigid/ fixed  sphere  of  radius  a 
centered  at  the  origin  and  an  incident  plane  wave  pmc  =  p0e~'kz  traveling  in  the  negative  z 

direction.  The  scattered  pressure  from  a  fluid  sphere  is  can  be  written  as  a  series  of  modal 
terms: 


PscaM’O’P)  =  TjBMcos  0)hH\kr), 


(2.17) 


where  Pm  is  the  Legendre  polynomial  and  /V '  ’  is  the  spherical  Hankel  function  of  the  first  kind 
(Anderson,  1950).  A  rigid/fixed  sphere  may  be  thought  of  as  a  fluid  sphere  with  infinite 
density.  In  this  case,  the  coefficient  Bm  becomes 


-p0HT(2rn  +  l)jm  m 
hfika) 


(2.18) 


where  jm  is  the  spherical  Bessel  function  of  the  first  kind.  The  derivatives  marked  by  the 
primes  are  with  respect  to  the  functions'  argument,  ka  .  For  backscattering,  6  -  0  and 
Pm(c os  6)  =  1  for  all  values  of  m  .  Simplifying,  the  backscattered  pressure  from  a  rigid/fixed 


sphere  is 
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(2.19) 


PiJn  = 


-Pq(-0,,'(2»7  +  \)hjkr)j,'(ka) 
™= o  h”’ ( ka ) 


The  backscattering  amplitude  and  target  strength  can  now  be  found  using  equations  (2.10)  and 

(2.11). 


2.3.3  Comparison  of  Kirchhoff  Method  and  Modal  Series  Solution  for  a 
Rigid/Fixed  Sphere 
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Figure  3  -  Target  strength  of  a  rigid/fixed  sphere  of  radius  5  cm,  from  the  Kirchhoff  method  (starred 
line)  and  modal  series  solution  (smooth  line). 


Target  strength  of  a  rigid/fixed  sphere: 
Kirchhoff  method  and  modal  series  solution 


4T  \ 


h —  Kirchhoff  integral 
- Modal  series  solution 


Figure  3  shows  how  the  Kirchhoff  method  result  seems  to  converge  to  the  exact  solution 
as  ka  increases.  However,  both  curves  have  strikingly  different  phases  for  their  peak  and  null 
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patterns.  The  reason  for  this  is  that  the  Kirchhoff  method  does  not  include  the  effects  of 
diffraction  around  the  sphere.  The  modal  series  solution  includes  both  the  effects  of  reflected 
and  diffracted,  or  creeping,  waves.  Described  by  Franz  (1957),  creeping  waves  diffract  around 
the  object  and  encircle  it  multiple  times.  Watson  (1918)  and  Sommerfeld  (1949)  described  the 
phenomenon  in  the  study  of  electric  waves  diffracting  around  the  Earth.  Creeping  waves 
constantly  radiate  off  energy  tangentially  from  the  object;  some  of  that  radiates  in  the 
backwards  direction  and  interferes  with  the  reflected  wave.  They  can  be  measured 
experimentally  by  using  short  sound  pulses  due  to  their  increased  travel  distance  and  time 
(Uberall  et  al.,  1966).  These  waves  can  also  be  influenced  by  surface  waves  if  the  object  is  elastic. 
With  a  rigid/ fixed  boundary  condition,  however,  they  only  depend  on  the  geometry  of  the 
scatterer.  As  ka  increases  and  these  effects  of  diffraction  lessen,  the  Kirchhoff  method  provides 
a  better  approximation  to  the  exact  modal  series  solution. 


2.4  Scattering  from  a  Finite  Cylinder 

In  this  section,  the  Kirchhoff  integral  is  solved  for  a  rigid / fixed  finite  cylinder.  The  results 
are  compared  to  an  approximate  modal  series  based  solution  for  a  rigid/ fixed  finite  cylinder. 
This  latter  approximation  is  based  on  the  exact  modal  series  solution  for  scattering  from  an 
infinite  cylinder.  The  approximation  neglects  the  effects  of  the  end  of  the  finite  cylinder,  and  is 
therefore  most  accurate  near  broadside.  Therefore,  the  Kirchhoff  integral  for  the  rigid/fixed 
finite  cylinder  is  solved  so  the  incident  wave  is  normal  to  the  cylinder's  axis. 
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2.4.1  Kirchhoff  Integral  Solution  for  a  Finite  Cylinder 


The  first  step  involves  defining  a  cylindrical  coordinate  system  with  radial  variable  r , 
polar  variable  0 ,  vertical  variable  z ,  and  the  same  Cartesian  coordinate  system  as  was  defined 
in  the  previous  section.  Assume  the  finite  cylinder  with  length  L  and  radius  a  is  centered  at 
the  origin  so  that  it  is  lengthwise  along  the  z  axis.  An  incident  plane  pinc  =  p0e,kx  wave 
traveling  in  the  positive  x  direction  strikes  the  cylinder  normal  to  its  axis.  Because  the  incoming 
wave  arrives  at  broadside,  the  cylinder  ends  are  not  insonified  and  therefore  do  not  contribute 
to  the  scattering  as  calculated  using  the  Kirchhoff  integral.  Therefore  only  the  sides  of  the 
cylinder  are  considered. 


Figure  4  -  Cylindrical  coordinate  system,  shown  with  incoming  wave  and  without  scatterer  geometry. 


It  can  thus  be  seen  that 

k=kx  (2.20) 

n  =  cos  Ox  +  sin  Oy 
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r'  =  a  cos  Ox  +  a  sin  Oy  +  zOz 
do  =  ad Odz . 

Substituting  these  expression  into  the  Kirchhoff  integral  gives 


.  LI 2  3^/2 
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The  integral  over  z  is  straightforward,  leaving 
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It  is  split  in  order  to  help  solve 
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Two  substitutions  are  made:  in  the  first  integral,  a  =  —0  +  n ;  in  the  second  integral,  fd  -  0  -  n  \ 
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Both  integrals  can  be  solved  with  help  of  the  relationship  (Gaunaurd,  1985) 


TC/l 


n , 


[  cos xe~izcosxdx  =  l--[tfj(z)  +  iJ,(z)\, 
i  2 


(2.25) 


where  H l  is  the  Struve  function  of  order  one,  and  is  the  Bessel  function  of  the  first  kind  of 
order  one.  Thus,  equation  (2.24)  becomes 

/  =  ^^[2  - n[Hx(2ka)  +  iJ{(2ka)]] .  (2‘26) 

2  n 


Target  strength  easily  follows  from  equation  (2.11). 
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2.4.2  Modal  Series  Based  Solution  for  a  Finite  Cylinder 


The  modal  series  approximation  for  the  scattering  of  finite  cylinders  is  based  on  the 
exact  modal  series  solution  for  an  infinite  cylinder.  It  assumes  that  end  effects  are  negligible, 
and,  therefore,  it  is  much  more  accurate  near  broadside  than  at  high  angles  of  incidence.  Stanton 
(1988)  gives  the  expression  for  scattering  from  a  rigid / fixed  finite  cylinder  of  length  L  and 
radius  a : 


f=  —  S^fjBm(-i)mcos(m<p) 


(2.27) 


n  A 


m= 0 


where 


A  = 


kL{r-r.)-rc 


(2.28) 


and 
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m 

nfm 


(2.29) 


The  r  terms  represent  the  unit  position  vectors  of  the  source  (subscript  i),  receiver  (subscript  r), 
and  cylinder's  axis  (subscript  c).  In  the  case  of  broadside  incidence,  A  =  0  .  The  term  <p  is  the 
angle  between  the  source  and  receiver  vectors  in  the  plane  perpendicular  to  the  cylinder  axis. 
For  the  backscattering  case,  (p  =  7T ,  and  therefore  cos (nup)  =  (-1)"' .  The  term  £m  is  equal  to  1  for 
m  =  0  and  equal  to  2  for  all  other  vales  of  m  .  Thus,  the  backscattering  amplitude  becomes 


iLy£m{-lTJm{ka) 
X  -0  Hfika) 


(2.30) 


Target  strength  follows  from  equation  (2.11). 
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2.4.3  Comparison  of  Kirchhoff  Method  and  Modal  Series  Based  Solution 


for  a  Rigid/Fixed  Finite  Cylinder 

Broadside  target  strength  of  a  rigid/fixed  finite  cylinder: 
Kirchhoff  method  and  modal  series  approximation 

-15  i - i - i - 


0  5  10  15 

ka 

Figure  5  -  Broadside  target  strength  of  a  rigid/fixed  finite  cylinder  of  length  10  cm  and  radius  1  cm 
using  the  Kirchhoff  method  (starred  line)  and  modal  series  approximation  (solid  line). 

As  with  the  sphere,  there  is  a  discrepancy  in  the  null-peak  phase  patterns  caused  by  the 
creeping  waves  not  accounted  for  in  the  Kirchhoff  method.  The  effects  of  this  are  again  lessened 
as  ka  increases.  There  seems  to  be  much  less  of  an  effect  from  creeping  waves  for  the  finite 
cylinder  than  for  the  sphere,  as  evidenced  by  the  former's  better  agreement  with  the  Kirchhoff 
method.  It  is  hypothesized  that  while  the  waves  diffract  relatively  well  around  the  entire  body 
of  the  sphere  and  the  curved  sides  of  the  cylinder,  they  do  not  diffract  as  well  around  the  flat 
ends  of  the  cylinder. 
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2.5  Kirchhoff  Method  Approximations  for  the  Sand  Dollar 


When  an  object's  geometry  becomes  too  complex,  it  may  not  be  possible  to  obtain  a 
closed  form  solution  to  the  Kirchhoff  integral.  However,  it  is  possible  under  some  circumstances 
to  approximate  the  complex  geometry  of  the  object  with  a  simpler  geometry  that  does  have  an 
analytic  solution  to  the  Kirchhoff  integral.  A  sand  dollar  is  an  example  of  a  complex  geometry 
where  the  integral  expression  cannot  be  solved  analytically  due  to  its  irregular  shape  and 
surface  roughness.  Generally,  sand  dollars  have  flat  bottoms  and  slightly  domed  tops  in 
addition  to  fine  scale  surface  roughness.  In  this  section,  the  Kirchhoff  integral  is  solved  for 
objects  that  approximate  both  of  the  orientations  of  the  sand  dollar.  A  round  disk  is  used  to 
model  the  flat  bottom  and  a  spherical  cap  is  used  to  model  the  rounded  top.  This  section  starts 
with  the  use  of  the  Kirchhoff  integral  to  derive  scattering  from  an  infinitely  thin,  rigid/ fixed 
disk,  as  this  sheds  light  on  the  derivation  for  the  (finite  thickness)  disk  and  spherical  cap. 


2.5.1  Kirchhoff  Integral  Solution  for  an  Infinitely  Thin  Disk 

The  first  step  in  solving  the  Kirchhoff  integral  for  a  disk  is  to  start  with  the  simplest  case, 
a  circular  rigid/  fixed  disk  of  radius  a  and  zero  thickness.  The  lack  of  thickness  means  that  only 
the  front  round  face  of  the  disk  is  incorporated  into  the  solution.  Returning  to  the  cylindrical 
coordinate  system,  assume  the  disk  is  located  at  the  origin,  perpendicular  to  the  z  axis.  The 

incoming  wave  vector  k  forms  an  angle  /?  with  the  z  axis,  varying  over  the  range  of  angles 
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Figure  6  -  Geometry  of  the  infinitely  thin  disk. 


In  this  coordinate  system,  the  relevant  parameters  of  the  Kirchhoff  integral  are  given  by 

&  =-£sin/?x  +  fccos/?z  (2.31) 

n  =  -z 

r  =  r  cos  Ox  +  r  sin  By 
do  -  rdrdO . 

Substituting  these  into  the  Kirchhoff  integral,  the  backscattering  amplitude  is  given  by 

-2i  cos  }  J  ^-mrsin/icose  -  -  -  (2>32) 


f  = 
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\r\e 
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dOdr . 


Abramowitz  and  Stegun  (1965)  give  an  integral  definition  of  a  Bessel  function  of  the  first  kind: 

(2.33) 


•n  ft 


J(z)  =  -\eizcos»cos(nO)dO. 
71  J 


Substituting  this  expression  back  into  equation  (2.32),  the  backscattering  amplitude  is 
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(2.34) 


/  =  lmcosP  j  rJ0(-2kr sin  J3)dr . 

A  o 

Next,  the  substitution  z  =  -2fcrsin  [5  is  made 

~  .  n  -2ka  sin  B  ,  .  n  ~2ka  sin  B 

,  -2m  cos  p  r  -z  ,  ,  .  -dz  -mcosp  r  T  .  .  , 

/= - 3 -  J  0/  •  •  R  =  0  3?  2  •  2  o  I  ^o(z)* 

/t  *  2x  sin  p  2k  sin  p  2Ak  sin  p  J0 

Abramowitz  and  Stegun  give  the  derivative  of  a  Bessel  function  as 

iiil 

Setting//  =  v  =  1 ,  this  simplifies  to 


(2.35) 


(2.36) 


~[zJl(z)]  =  zJ0(z). 
dz 


(2.37) 


Substituting  and  solving. 


-mcosp  r  ,  iaJA-2ka  sin p) cos p 

/  =  o3?2  •  2  o  I  ^L^i(z)J  = - 1 - ^—r, - 

2Ak  sin  P  J()  2  sin  p 


(2.38) 


The  target  strength  of  the  rigid/ fixed,  infinitely  thin  disk  is  thus  given  by 


TS  = 10 log 


aJ  j  ( — 2ka  sin  P)  cos  P 


2  sin  P 


(2.39) 


Urick  (1983)  gives  this  result  in  a  slightly  different  form  as  the  target  strength  of  a  circular  plate. 


2.5.2  Kirchhoff  Integral  Solution  for  a  Disk  with  Finite  Thickness 

The  next  step  in  solving  the  Kirchhoff  integral  for  the  disk  is  to  give  it  a  thickness.  There 
are  now  two  different  scattering  regions  of  the  disk:  the  flat  round  face  and  the  curved  edge. 
The  height  of  this  new  edge  is  the  disk's  thickness,  T. 
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Figure  7  -  Geometry  of  the  disk. 

The  disk  is  simply  a  finite  cylinder  with  an  extremely  high  radius  to  length  aspect  ratio. 
Therefore,  this  derivation  is  like  that  of  the  finite  cylinder  over  any  angle  of  incidence.  The 
Kirchhoff  integral  can  be  split  to  analyze  each  of  the  disk's  elements.  These  two  regions  of 
scattering  add  coherently  to  give  the  total  scattering  amplitude, 

/=/,«.+/*■  <Z4°) 

The  backscattering  amplitude  of  the  face  is  very  similar  to  before,  with  the  only  difference  being 

a  phase  shift  due  to  the  face's  negative  displacement  by  T/2, 

_ika2Jx(-2ka  sin  f3)  cos /3  -ikTcosp  (2-41) 

face  2ka  sin  j3 

Next,  the  scattering  from  the  edge  is  evaluated.  The  relevant  parameters  within  the  Kirchhoff 
integral  are  given  by 

k  =  —k  sin  /3.x  +  k  cos  /3z  (2.42) 
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n  =  cos  Ox  +  sin  Oy 


r'  =  a  cos  Ox  +  a  sin  Oy  +  zz 
cl  <7  =  ad  Odz . 

Substituting  into  the  integral  expression,  the  backscattering  amplitude  is  given  by 


.  772  71 1 2 


f edge  ~  J"  J  -Sin^COS^C 


il(—ka  sin  (5  cos  O+kz  cos  J3 ) 


ad  Odz  ■ 


(2.43) 


-772 -/r/2 


The  integral  over  z  is  straightforward: 


/« 


-ia  sin(^r  cos  J3)  sin  /3 


k!  2 


[  cos  Oe  1 

-nil 


sin  ficosd 


dO. 


(2.44) 


2n  cos  j3 

To  solve  this  integral  it  first  must  be  split.  Then  a  substitution  a  =  —0  is  made  in  one  of  the 
resulting  integrals: 


f 

J  edge 


-ia  sin (kT  cos  /?)  sin  / 3 


f edge 


27TCOS 

-ia  sin  (kT  cos  j3)  sin  ft 


2^cos  P 


0 

nil 

f  cos  Oe->2kasARpmsedO  + 

J  cos  Oe'i2kasinflcosed0 

-nil 

0 

nil 

nil 

J  cos  ae-i2kasiD/}cosada+ 

f _ /x  -ilkasin  BcosO  i 

I  cosOe  H  ad 

0 

0 

(2.45) 


(2.46) 


These  integrals  may  now  be  solved  with  help  from  equation  (2.25),  giving 
-ia  sin  (kT  cos  p)  sin  P  i 


f 

J  edge 


-  2  -n\Hx(2ka  sin p)  +  iJ x{2ka  sin /?)]]  . 


(2.47) 


2^cos  P 

Gaunaurd  (1985)  presents  a  similar  derivation  and  arrives  at  the  same  result  in  his  discussion  of 
total  insonification  of  a  finite  cylinder.  The  target  strength  of  the  rigid/ fixed  disk  is  thus  given 
by 
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re  =  io log  i  f.  +  f,  i 


=  lOlog 


face  J  edge 

ika2Jf—2ka  sin /?) cos /3  ^-ikTcosp  ia  sin(&r cos /?)sin  /3 


(2.48) 


2ka  sin  [3  2k  cos  [3 

x  [2  -  K  [//,  (2ka  sin  J3)  +  ij  f2ka  sin  /?)]] 


2.5.3  Kirchhoff  Integral  Solution  for  a  Spherical  Cap 

The  next  object  that  the  Kirchhoff  integral  is  solved  for  is  a  spherical  cap,  modeling  the 
sand  dollar's  top  side.  Returning  to  the  spherical  coordinate  system,  assume  a  sphere  of  radius 
a  centered  at  the  origin.  A  slice  perpendicular  to  the  x-y  plane  is  made  through  the  top  half  of 
the  sphere,  and  everything  below  is  discarded,  leaving  only  the  top  cap  of  the  sphere.  Another 
way  to  define  the  cap  is  by  polar  angle,  0  =  0  to  some  angle  rj .  The  incident  wave  makes  an 
angle  [3  with  the  z  axis. 


z 


Figure  8  -  Geometry  of  the  spherical  cap. 


The  relevant  parameters  of  the  Kirchhoff  integral  are  given  by 
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k  =  —k  sin  fix  +  k  cos  fiz 


(2.49) 


n  =  cos  (p  sin  Ox  +  sin  <p sin  Oy  +  cos  Oz 
r'  =  a  cos  ^sin  Ox  +  a  sin  ^sin  Oy  +  a  cos  Oz 
da  =  a1  sin  Od  Od  (p . 

Substituting  into  the  Kirchhoff  integral  and  simplifying, 

1  In  11 


f  f  [  (sin /?  cos  ^  sin# -cos/?  cos#)  sin  #e 

O  77*  J  J 


i  2  ka  (sin  fi  cos  (p  sin  9— cos/?  cos 


2  71 


e)d  Od  (p . 


(2.50) 


o  o 


This  author  is  not  aware  of  an  analytic  solution,  so  this  double  integral  is  to  be  evaluated 
numerically. 


2.6  Solving  the  Kirchhoff  Integral  Numerically 

Using  a  simplified  approximation  to  the  geometry  of  a  complex  shape  in  order  to  obtain  a 
closed  form,  analytic  solution  to  the  Kirchhoff  integral  can  introduce  errors  into  the  predicted 
scattering  associated  to  the  geometric  differences.  In  cases  where  analytic  simplification  is 
impractical,  or  the  resulting  errors  are  unacceptable,  it  is  possible  to  use  a  numerical  approach 
to  solve  the  Kirchhoff  integral  that  does  not  rely  on  simplifying  the  object's  shape.  In  this 
section,  a  method  is  presented  for  solving  the  Kirchhoff  integral  numerically,  regardless  of  the 
object's  irregularities. 
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2.6.1  Creating  Geometry  and  a  Mesh  in  COMSOL  Multiphysics 

The  first  step  in  solving  the  Kirchhoff  integral  numerically  is  to  create  a  surface  mesh  of 
the  object  of  interest.  This  is  simply  a  collection  of  small  connected  triangles  that  attempt  to 
match  the  surface  of  the  object.  The  size  of  these  triangles  determines  how  well  they  conform  to 
the  shape.  One  simple  method  for  creating  the  object's  surface  mesh  is  through  use  of  the  finite 
element  software  COMSOL  Multiphysics®  (ver.  3.2a),  referred  to  in  this  work  simply  as 
COMSOL.  COMSOL  is  a  commercially  available  program  for  solving  finite  element  method 
problems,  which  will  be  explored  further  in  Chapter  3.  For  now,  it  is  used  for  its  capability  to 
model  and  mesh  complex  geometries. 

Within  COMSOL,  the  geometry  of  the  object  of  interest  can  be  created  using  either  the 
GUI  or  script.  Objects  can  be  created  from  scratch  or  imported  as  computer-aided  design  (CAD) 
files.  Once  the  geometry  is  created,  a  mesh  can  be  generated  to  conform  to  it.  The  mesh  for  a 
three  dimensional  object  consists  of  tetrahedrons  in  the  interior  and  triangles  on  the  surfaces. 
The  vertices  of  these  shapes  are  known  as  mesh  vertices.  The  mesh  is  created  by  defining  input 
parameters  that  regulate  the  size  of  the  triangles  and  tetrahedrons  and  how  well  they  conform 
to  curved  parts  of  the  geometry.  The  size  of  the  triangles  on  the  geometry's  surface  will  have  an 
effect  on  the  accuracy  of  numerical  solution  to  the  integral,  since  each  triangle  acts  a  differential 
element. 

The  mesh  may  be  represented  by  two  sets  of  information:  coordinates  and  connectivity. 
The  coordinate  information  lists  the  location  of  each  of  the  numbered  mesh  vertices.  The 
connectivity  lists  the  vertices  (by  number)  that  make  up  each  of  the  mesh  elements,  in  some 
systematic  way  (e.g.  counterclockwise  for  triangles,  right  hand  rule  for  tetrahedrons).  In  this 
way,  each  mesh  tetrahedron  is  defined  by  a  connectivity  of  four  vertices  whose  locations  are 
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contained  in  the  coordinate  list.  Similarly,  each  triangle  is  defined  by  three  vertices.  Because  the 
Kirchhoff  method  uses  a  surface  integral,  only  the  information  regarding  the  surface  triangles  is 
important.  This  can  easily  be  separated  and  assembled  within  COMSOL  into  two  arrays. 


2.6.2  Determination  of  Insonified  Surface 

Once  the  two  surface  mesh  arrays  have  been  compiled,  the  next  step  is  to  determine 
which  triangular  surface  elements  should  contribute  to  the  solution.  Because  the  Kirchhoff 
integral  is  only  applied  to  the  insonified  region  of  the  object,  any  surface  in  the  shadow  needs  to 
be  ignored.  When  solving  the  integral  analytically,  the  proper  surface  of  the  object  is  defined  by 
the  limits  of  integration.  Since  these  are  not  defined  here,  there  needs  to  be  some  way  to 
eliminate  the  triangle  elements  that  are  not  in  the  line  of  sight  of  the  incoming  wave  vector.  This 
problem  is  similar  to  one  in  computer  graphics,  where  several  surfaces  or  polygons  share  the 
same  set  of  pixels.  In  order  to  display  the  correct  surface  (e.g.  that  which  is  closest  to  the 
viewer),  a  process  known  as  hidden  surface  removal  determines  what  information  each  pixel 
should  display  (De  Berg  et  ah,  1997).  Newell's  algorithm  is  a  similar  method  for  determining 
which  polygons  are  hidden  (Newell,  1972).  A  simplified  version  of  Newell's  algorithm  is 
developed  here  in  order  to  determine  which  triangles  are  hidden  by  others,  from  the  view-point 
of  the  incoming  wave  vector. 

To  use  the  method  developed  here,  a  step  must  be  taken  before  the  mesh  is  created.  In 
order  to  avoid  complicated  coordinate  transforms  later  on,  the  incoming  wave  vector  is  set  to 
always  coincide  with  the  negative  z  direction.  The  object  can  be  rotated  in  COMSOL  so  that  the 
desired  angle  of  incidence  is  achieved.  After  the  mesh  is  created  for  this  orientation  of  the 
object,  the  simplified  version  of  Newell's  algorithm  can  be  applied: 
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1)  Sort  elements  by  max  z  value 

2)  Minmax  test  in  x 

3)  Minmax  test  in  y 

4)  Rasterization  test 

For  the  first  step,  the  elements  can  be  sorted  by  their  maximum  z  values  by  comparing  the 
positions  of  each  triangle's  three  vertices.  Because  the  incoming  wave  vector  travels  in  the 
negative  z  direction,  the  elements  higher  on  this  list  are  more  likely  to  be  part  of  the  front 
surface. 

The  remainder  of  the  algorithm  involves  determining  whether  elements  are  covered  by 
others.  If  the  projections  into  the  x-y  plane  of  two  elements  intersect,  then  the  element  with  the 
greater  z  value  covers  the  other.  Elements  which  are  covered  by  others  must  lie  on  a  surface  that 
is  not  in  the  line  of  sight  of  the  incoming  wave  vector  and  therefore  in  the  shadow.  These 
hidden  elements  are  removed  from  the  list  as  they  are  found.  The  tests  proceed  by  starting  with 
the  first  element  in  the  list  and  comparing  it  to  all  the  other  elements  below  it  in  the  list.  When 
these  comparisons  are  finished  and  covered  elements  are  removed,  the  second  element  in  the 
list  is  compared  to  all  the  other  elements  below  it.  This  process  eliminates  redundant 
comparisons. 

The  first  test  is  a  minmax  test  in  x:  determining  whether  the  two  elements  share  any  x 
values  by  comparing  minimum  and  maximum  values.  If  this  test  fails  (they  do  not  share 
common  x  values),  the  second  element  is  not  covered  by  the  first,  and  the  next  element  is 
compared. 
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Figure  9  -  Examples  of  minmax  tests  in  x  for  Newell's  algorithm. 


If  the  two  elements  pass  the  minmax  tests  in  x  and  y  and  share  values  in  both  x  and  y,  a 
rasterization  test  is  performed.  Given  a  sampling  rate,  a  rectangle  of  x-y  points  is  created;  the 
coordinates  of  these  points  are  determined  by  the  minimum  and  maximum  x  and  y  values  of 
the  first  element.  Each  point  is  sampled  to  determine  if  it  is  part  of  either  triangular  element. 
Points  on  the  border  between  two  triangles  are  considered  part  of  neither.  If  a  point  is  part  of 
both,  then  the  first  element  must  cover  the  second,  which  is  then  removed  from  the  list.  This  test 
is  run  only  after  a  set  of  elements  has  passed  both  minmax  tests  to  save  run  time. 


2.6.3  Evaluating  the  Kirchhoff  Integral  Numerically 

With  the  final  list  of  elements  representing  the  object's  front  surface,  the  Kirchhoff 
integral  can  now  be  evaluated.  There  is  some  ambiguity  in  the  direction  of  the  normal  vectors 
that  must  be  clarified.  Due  to  the  numbering  scheme  in  the  COMSOL  connectivity  list. 
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determining  the  cross  product  of  the  triangle  edges  does  not  always  provide  an  outward 
normal.  This  can  be  fixed  by  noting  that  the  dot  product  of  the  normal  and  unit  wave  vectors 
should  either  be  negative  or  zero,  because  the  faces  are  either  oriented  toward  the  positive  z 
direction  or  are  perpendicular  to  it.  None  of  the  triangles  can  be  facing  the  negative  z  direction 
after  applying  Newell's  algorithm.  For  the  elements  perpendicular  to  the  wave  vector,  it  is  not 
important  if  the  normal  vector  faces  inward  or  outward  since  the  dot  product  and  integrand 
reduces  to  zero  in  either  case.  For  the  other  faces  however,  it  is  simply  enough  to  compute  the 
absolute  value  of  the  dot  product  of  the  normal  and  unit  wave  vectors  and  set  it  negative.  For 

the  kinc  ■  r'  term  in  the  exponent,  the  average  of  the  dot  products  of  the  wave  vector  and 

position  vectors  of  each  vertex  is  used.  The  area  of  each  triangle  is  computed  for  the  differential 
element  term.  Finally,  the  expression  is  summed  over  each  triangle  element  and  the  Kirchhoff 
integral  can  be  numerically  evaluated. 


2.7  Comparison  of  Methods  for  Solving  the  Kirchhoff 
Integral 

The  numerical  method  of  solving  the  Kirchhoff  integral  can  be  tested  against  some  of  the 
analytic  solutions  derived  earlier.  In  each  case,  the  mesh  is  made  fine  enough  so  that  there  is 
good  agreement  between  the  two  methods.  The  mesh  size  at  which  the  numerical  solution 
converged  to  the  analytic  solution  varied  in  each  case,  and  proved  highly  dependent  on  the 
geometry  of  the  scatterer.  Therefore,  a  simple  rule  for  determining  an  appropriate  mesh  size 
was  not  ascertained. 
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2.7.1  Sphere 


A  sphere  of  radius  5  cm  was  created  in  COMSOL  and  meshed,  seen  in  figure  10.  Because 
of  its  symmetry,  no  care  needed  to  be  taken  to  orient  the  sphere  in  any  fashion. 


Figure  10  -  Surface  mesh  of  a  sphere  of  radius  5  cm  in  COMSOL. 


This  mesh  had  a  maximum  element  size  of  0.5  cm,  giving  4266  surface  triangle  elements. 
Extracting  the  connectivity  and  node  coordinate  information,  these  elements  were  subjected  to 
the  Newell's  algorithm  in  the  technical  computing  software  MATLAB®  (ver.  7.2.0.232).  To 
visualize  the  effect  of  this,  the  remaining  triangles  are  seen  plotted  in  figure  11. 
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Figure  11  -  Surface  mesh  of  a  sphere  in  MATLAB  after  Newell's  algorithm. 

Newell's  algorithm  has  the  desired  effect  of  removing  the  elements  out  of  the  line  of  site  of  the 
incoming  wave  vector  from  above.  To  ensure  that  each  of  the  triangles  in  the  plot  were 
represented  by  an  element  and  not  simply  formed  by  surrounding  triangles,  the  midpoint  of 
each  were  plotted,  seen  in  figure  12. 
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Figure  12  -  Surface  mesh  of  a  sphere  in  MATLAB  after  Newell's  algorithm,  with  triangle  midpoints. 


Taking  a  closer  look,  it  is  clear  that  Newell's  algorithm  kept  all  the  desired  triangular  elements. 


Figure  13  -  Close  view  of  surface  mesh  of  sphere,  with  midpoints. 
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Finally,  the  Kirchhoff  integral  was  solved  numerically  over  each  element  and  compared  to  the 
analytic  solution  for  a  rigid/ fixed  sphere  (section  2.3.1)  for  different  values  of  ka,  the  product  of 
the  wavenumber  and  the  sphere's  radius,  shown  in  figure  14. 
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Figure  14  -  Target  strength  of  a  rigid/fixed  sphere  with  radius  5  cm  by  solving  the  Kirchhoff  integral 
analytically  (solid  line)  and  numerically  (stars). 


By  using  a  sufficiently  fine  surface  mesh,  the  numerical  solution  converged  to  the  analytic 
solution. 
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2.7.2  Finite  Cylinder 


A  finite  cylinder  of  length  10  cm  and  radius  1  cm  centered  at  the  origin  was  created  in 
COMSOL,  oriented  with  its  lengthwise  axis  along  the  x  axis  so  that  the  incoming  wave  vector, 
along  the  negative  z  axis,  met  it  at  broadside.  The  cylinder  was  then  meshed  with  a  maximum 
element  size  of  0.15  cm,  producing  11170  surface  elements. 


Figure  15  -  Surface  mesh  of  a  finite  cylinder  in  COMSOL. 

Again,  the  triangular  elements  were  run  through  Newell's  algorithm  in  MATLAB,  with  the 
remaining  ones  shown  in  figure  16. 
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Figure  16  -  Surface  mesh  of  a  finite  cylinder  in  MATLAB  after  Newell's  algorithm. 

Note  that  the  elements  on  the  ends  of  the  cylinder  were  not  removed.  This  is  due  to  a 
shortcoming  in  the  rasterization  test  in  Newell's  algorithm.  Because  the  elements  are  vertical, 
they  cannot  be  sampled  in  the  x-y  plane.  However,  because  the  dot  product  of  the  unit  normal 
and  unit  wave  vector  is  zero,  their  contribution  amounts  to  zero  regardless.  The  Kirchhoff 
integral  can  now  be  evaluated  over  all  the  remaining  elements  and  compared  to  the  analytic 
solution  for  a  finite  cylinder  (section  2.4.1)  over  a  range  of  values  of  ka,  the  product  of  the 
wavenumber  and  the  cylinder's  radius.  Again,  the  numerical  results  converge  to  the  analytic 
solution,  as  seen  in  figure  17. 
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Comparison  of  methods  for  solving  Kirchhoff  integral 
for  a  finite  cylinder,  broadside  incidence 
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Figure  17  -  Target  strength  of  a  rigid/fixed  finite  cylinder  with  length  10  cm  and  radius  1  cm  by  solving 
the  Kirchhoff  integral  analytically  (solid  line)  and  numerically  (stars). 


2.7.3  Disk 


A  disk  of  radius  5  cm  and  thickness  1  cm,  centered  at  the  origin,  was  created  in  COMSOL, 
oriented  so  that  the  incoming  wave  vector  along  the  negative  z  axis  arrived  normal  to  the  round 
face.  The  surface  mesh  had  a  maximum  element  size  0.25  cm,  giving  10640  surface  triangles. 
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Figure  18  -  Surface  mesh  of  a  disk  in  COMSOL. 


The  triangular  elements  were  then  run  through  Newell's  algorithm  in  MATLAB,  with  the 
results  shown  in  figure  19. 


Figure  19  -  Surface  mesh  of  a  disk  in  MATLAB  after  Newell's  algorithm. 
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Many  of  the  vertical  elements  on  the  side  of  the  disk  remain,  as  they  did  on  the  ends  of  the  finite 
cylinder;  however,  they  do  contribute  to  the  solution.  Finally,  the  Kirchhoff  integral  can  be 
solved  numerically  over  all  the  remaining  elements,  and  the  solution  is  compared  to  the 
analytic  solution  (section  2.5.2).  Again  the  numerical  results  matched  the  analytic  solution,  as 
seen  in  figure  20. 
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Figure  20  -  Target  strength  of  a  rigid/fixed  disk  with  radius  5  cm  and  thickness  1  cm  by  solving  the 
Kirchhoff  integral  analytically  (solid  line)  and  numerically  (stars). 
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2.7.4  Spherical  Cap 


To  generate  the  spherical  cap,  a  sphere  of  radius  10  cm  was  first  created  in  COMSOL, 
centered  at  the  origin.  The  top  30  degrees  of  the  sphere  were  kept,  leaving  a  cap  of  radius  5  cm 
and  thickness  1.34  cm.  The  cap  was  oriented  so  that  the  incoming  wave  along  the  negative  z 
axis  arrived  normal  to  its  top.  It  was  meshed  with  a  maximum  element  size  of  0.25  cm,  giving 
9382  surface  elements. 


Figure  21  -  Surface  mesh  of  a  spherical  cap  in  COMSOL. 

The  triangular  elements  were  then  run  through  Newell's  algorithm  in  MATLAB,  with  the 
results  plotted  in  figure  22. 
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Figure  22  -  Surface  mesh  of  a  spherical  cap  in  MATLAB  after  Newell's  algorithm. 

Finally,  the  Kirchhoff  integral  was  solved  numerically  over  all  the  remaining  elements  and 
compared  to  the  analytic  solution  for  a  spherical  cap  (section  2.5.3)  over  a  range  of  ka,  where  a 
is  the  radius  of  the  cap.  Again,  the  numerical  results  match  the  analytic  solution,  as  seen  in 
figure  23. 
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Comparison  of  methods  for  solving  Kirchhoff  integral 
for  a  spherical  cap,  normal  incidence 


Figure  23  -  Target  strength  of  a  rigid/fixed  spherical  cap  with  radius  5  cm  and  radius  1.34  cm  by 
solving  the  Kirchhoff  integral  analytically  (solid  line)  and  numerically  (stars). 


2.8  Chapter  Summary 


In  this  chapter,  the  Kirchhoff  method  has  been  explored  for  understanding  scattering 
from  objects  with  and  without  known  exact  solutions.  The  derived  solutions  from  this  method 
have  their  limitations  and  may  not  be  valid  in  many  situations,  depending  on  factors  such  as 
boundary  condition  and  incident  wavelength.  However,  as  will  be  seen,  this  method  proves  to 
be  useful  in  studying  the  scattering  from  sand  dollars  at  certain  frequencies  and  angles  of 
orientation.  The  analytic  solutions  for  a  disk  and  spherical  cap  are  compared  to  experimental 
scattering  data  from  a  sand  dollar  in  Chapter  4.  The  numerical  approach  to  solving  the 
Kirchhoff  integral  was  presented  in  cases  where  representing  a  complex  geometry  with  a 
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simpler  one  introduces  too  many  errors.  This  approach  has  been  shown  to  agree  very  well  with 
analytic  solutions  to  the  Kirchhoff  integral.  Also  in  Chapter  4,  the  scattering  from  a  rigid/fixed 
three  dimensional  model  of  the  sand  dollar  is  calculated  with  this  numerical  method  and 
compared  to  experimental  results. 
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Chapter  3  The  Finite  Element  Method 


3.1  Introduction 

In  this  chapter,  the  theory  behind  the  finite  element  method  and  its  use  in  acoustic 
scattering  problems  are  explained.  This  is  presented  in  the  form  of  a  general  problem  of  acoustic 
scattering  from  a  rigid/ fixed  object.  A  technique  for  reducing  random  numerical  error  in  the 
solution  of  two  and  three  dimensional  acoustic  scattering  problems  is  then  presented.  COMSOL 
Multiphysics,  a  program  for  solving  finite  element  method  problems  is  explored.  A  particular 
focus  concerns  the  radiation  condition  used  by  this  program.  The  finite  element  method  is 
tested  for  scattering  from  a  rigid/fixed  infinite  cylinder  and  a  rigid/fixed  sphere,  problems  with 
known  exact  solutions.  The  methodology  for  modeling  the  sand  dollar  using  the  finite  element 
method  is  then  discussed. 


3.2  Theoretical  Background 

The  finite  element  method  is  a  general  numerical  technique  for  solving  boundary  value 
problems.  Its  application  to  acoustic  scattering  problems  is  now  shown  by  introducing  a  general 

three  dimensional  problem.  Consider  a  harmonic  plane  wave  pmc  =  p0e,k  x  that  is  incident  upon 

a  rigid/ fixed  object  in  free  space,  represented  by  Q  .  The  surface  of  the  object  is  labeled  T . 
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n 


Figure  24  -  Domain  of  three  dimensional  scattering  problem,  shown  with  incident  wave. 


The  governing  equation  for  the  total  pressure  field  in  the  domain  Q.  from  the  object  is  the 
Helmholtz  equation,  the  time  independent  form  of  the  wave  equation, 

V2p  +  k2p  =  OinH.  (3-1) 

The  total  pressure  is  the  sum  of  the  incident  and  scattered  pressure  fields, 

P  =  Pine +  P  scaf  (3‘2) 

Therefore,  the  scattered  pressure  field,  pscat ,  also  satisfies  the  Helmholtz  equation  in  , 

V2Psca,+klPsca,=  0  in  ^ 

The  goal  of  the  problem  is  to  determine  the  scattered  pressure  field. 

On  the  boundary  of  the  object,  T,  the  rigid/fixed  condition  requires  that  the  normal 
gradient  of  total  pressure  is  zero, 

n  ■  Vp  =  0,  (3.4) 


where  n  is  the  surface  normal  of  the  object.  This  is  easily  written  for  the  scattered  pressure  in 
terms  of  the  incident  pressure, 

«  •  ^ Pscat  =  -«  '  VP,nc  ■  (3-5) 
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To  complete  the  boundary  value  problems,  a  condition  of  the  behavior  of  the  scattered  pressure 
at  infinity  is  required,  given  by  the  Sommerfeld  radiation  condition  (Sommerfeld,  1949)  in  three 
dimensions  as 


dr 


scat 


-  ikp 


scat 


=  0. 


(3.6) 


This  condition  describes  the  behavior  of  the  scattered  pressure  as  r ,  the  distance  from  the 
scatterer,  goes  to  infinity. 

Note  that  the  problem  above  is  defined  on  the  unbounded  domain  Q  .  The  finite 
element  method  used  in  this  thesis  requires  a  domain  where  the  boundary  value  problem  is 
defined  to  be  bounded.  Therefore,  the  Sommerfeld  condition  is  substituted  with  an  approximate 
radiation  condition  defined  on  the  boundary  rinf  that  contains  the  object.  The  geometry  of  this 
new  bounded  region,  labeled  Q ,  is  seen  in  figure  25. 


Figure  25  -  Geometry  of  a  bounded  three  dimensional  scattering  problem. 


At  rinf/  the  condition  used  instead  of  equation  (3.6)  is  defined  as 
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(3.7) 


%2L=ap,..,  on  r, 

or 

and  equation  (3.3)  is  defined  over  Q  .  Therefore,  the  boundary  value  problem  that  will  be 
solved  with  the  finite  element  method  here  is  summarized  as 

V2Psca,  +k2Psca,  =  0  in  (3‘8) 

»  '  VPu:a,  =  ■  VPinc  °n 

%2L  =  *p,„  »nr„. 

or 

To  apply  the  finite  element  method,  the  boundary  value  problem,  given  in  (3.8),  is 
written  in  the  form  of  a  variation  equation,  as  shown  in  the  following  development.  First  the 
condition  on  Q  is  multiplied  by  a  weight  function  we  V  and  integrated  over  the  domain, 
resulting  in 

MvT,,,+*7>,„,>m=0  <3'9) 

n 

V  is  a  linear  space  of  functions  which  will  be  defined  in  a  moment,  and  the  asterisk  denotes  the 
operation  of  conjugation.  The  first  term  on  the  left  side  on  equation  (3.9)  is  integrated  by  parts, 
using  the  identity 

div(fVg)  =  Vf-Vg  +  fV2g.  (310) 

This  identity  is  valid  for  any  two  scalar  fields  /  and  g  which  are  twice  differentiable. 
Integrating  this  identity  over  the  domain  Q.  and  using  the  divergence  theorem  (Strang,  1991) 
gives 

J  fV2gdQ  =  - J  V/  ■  VgdQ  +jfVg ■  ddQ ,  (311) 

n  n  an 
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where  n  is  the  unit  normal  vector  outward  from  the  boundary  of  Q. .  Note  that  the  right  side  of 
this  expression  requires  /  and  g  to  be  only  first-order  differentiable.  Using  equation  (3.11)  in 

equation  (3.9)  with  /  substituted  with  w*  and  g  by  pscat  gives 


J (-Vw*  ■  Vpscat  +  k2w*pscat ) dQ  +  J  w\pscal  •  MY  =  0  Vw e  V  .  (312) 

n  an 

Using  dfl  =  r  +  rinf  and  the  boundary  conditions  in  (3.8)  gives  the  following  variational 


equation  for  the  scattered  field  pscat  e  V 


J  (Vvv  ■  Vpscat  -  k 2  w  pscat  )d£>-  J  ikw  psca,dY inf  =  -J  w\pinc  ■  MY  (313) 


Vwe  V . 

From  this  equation,  the  conditions  on  the  space  V  are  deduced:  the  functions  have  to  be  square- 
integrable  with  square-integrable  gradients.  Equation  (3.13)  can  be  written  in  compact  form  by 
introducing  the  notation 


“(W’Psca,)  =  j (Vw* '  Vpa 


Y(w)  =  -  J  w\pjnc  ■  MY 

r 


(3.14) 

(3.15) 


so  that  the  variational  equation  is 

a(yv,pKat)  =  y(w)  VweV. 


(3.16) 


Note  that  the  form  is  linear  with  respect  to  the  first  and  second  arguments  and  that  /(■)  is 
linear  with  respect  to  its  argument  as  well. 

The  finite  element  method  is  a  numerical  technique  to  construct  approximate  solutions 
to  variational  problems  such  as  (3.16).  An  approximation  to  the  infinite-dimensional  space  V  is 
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constructed.  This  space  is  denoted  Vh  and  has  the  property  that  it  is  a  subspace  of  V  .  The 
problem  is  solved  in  this  new  space,  that  is,  the  goal  is  to  determine  ph  such  that 

a(w",ph)  =  r(wh)  \/wh  e  Vh .  (3-17) 

The  space  V h  is  a  finite  dimensional  space,  which  means  that  a  basis  can  be  constructed 
using  a  finite  number  of  terms.  The  elements  of  that  basis  are  denoted  ftt(x) ,  with  i  e  { 1 , . . . ,  /V } . 

N  is  the  number  of  elements  in  the  basis  as  well  as  the  dimension  of  the  space  V h .  Therefore, 
ph  and  wb  may  be  represented  as 


Ph(x)  =  ^jpifti(x)/ 


1=1 


(3.18) 


N 

w‘(*)  =  j>^(x). 


(3.19) 


/= i 

Substituting  these  expressions  into  equation  (3.17)  and  using  the  linearity  of  and  /(■) 
results  in  the  algebraic  equation 


w-Kp  =  w-F  \/weRN. 


(3.20) 


The  term  p  =  (pl,...,pN)T  is  the  vector  components  of  ph  in  the  basis  ft ,  and  likewise  for  w . 
This  condition  is  satisfied  if  and  only  if  the  vector  p  is  the  solution  of  the  matrix  equation 

Kp  =  F .  (321) 


The  components  of  the  matrix  K  and  the  vector  F  are  defined  as 

Ky  =  aift^ftj),  (3.22) 

F  =  r(ft,) 

To  solve  equation  (3.21),  direct  or  iterative  solvers  are  used  (Golub  and  Van  Loan,  1996). 
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The  finite  element  method  discretizes  the  domain  of  a  problem  into  small  connected 
subdomains  of  some  polygonal  shape  (usually  triangles  or  tetrahedrons).  This  pattern  of 
discretization  is  referred  to  as  the  domain's  mesh.  The  individual  subdomains  within  the  mesh 
are  known  as  the  finite  elements  themselves.  Each  finite  element  is  formed  by  nodes  on  its 
boundary,  and  elements  are  related  by  their  shared  nodes.  If  an  element  only  has  nodes  at  its 
vertices,  it  is  said  to  be  linear.  Quadratic  elements  are  formed  by  designating  the  midpoints  of 
the  line  segment  between  the  vertices  as  nodes.  Even  higher  order  elements  may  be  formed 
with  additional  nodes,  but  their  utility  is  limited  with  the  techniques  available  (Becker  et  al, 


1981). 


Figure  26  -  Examples  of  two  linear  (left)  and  three  quadratic  (right)  triangular  finite  elements,  used  to 
discretize  two  dimensional  problems.  The  black  dots  represent  the  positions  of  the  nodes. 

In  this  work,  quadratic  triangular  Lagrange  elements  were  chosen  for  the  discretization  of  the 
two  dimensional  domains  and  quadratic  quadrilateral  Lagrange  elements  for  the  three 
dimensional  domains.  Quadratic  quadrilateral  elements  are  composed  of  ten  nodes,  one  for 
each  of  the  four  vertices  and  one  on  each  midpoint. 
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Figure  27  -  Example  of  a  three  dimensional  quadrilateral  quadratic  finite  element.  The  black  dots 
represent  the  positions  of  the  ten  nodes. 

The  basis  functions  are  chosen  based  on  the  nodes  of  the  problem.  In  this  work,  each 
basis  function  $  is  a  piecewise  quadratic  function,  equal  to  one  at  node  i  and  equal  to  zero  at 
all  other  nodes.  With  this  choice  of  basis,  the  solutions  for  scattered  pressure  at  the  nodes  are 
conveniently  contained  within  the  vector  p  . 

The  major  simplification  in  the  above  derivation  concerns  the  outer  boundary  on  which 
the  simplified  Sommerfeld  radiation  condition  was  assigned.  The  difficulty  in  solving  scattering 
problems  with  the  finite  element  lies  with  trying  to  reconcile  infinity  with  a  finite  numerical 
domain.  Because  the  goal  in  many  cases,  as  it  is  here,  is  to  understand  the  scattering  from  an 
object  located  in  free  space,  away  from  boundaries  and  obstacles,  scattering  problems  usually 
take  place  in  an  infinite  domain.  Since  it  is  impossible  to  create  an  infinite  numerical  domain  in 
the  finite  element  method,  such  a  problem  must  be  simplified  to  only  the  scatterer  and  a  finite 
region  surrounding  it.  The  scattered  pressure  can  be  determined  anywhere  within  this 
surrounding  region  through  the  finite  element  method,  but  not  outside  of  it.  A  technique  for 
getting  around  this  shortcoming  is  presented  in  the  next  section. 
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The  exterior  of  the  surrounding  medium  region  must  be  given  a  boundary  condition  that 
represents  free  space  extending  to  infinity,  known  as  a  radiation  condition.  This  boundary  must 
model  the  effects  of  the  free  space  it  represents;  specifically,  it  must  not  reflect  outgoing  waves 
back  into  the  domain.  Local  absorbing  boundary  conditions  (Atkinson,  1949;  Wilcox,  1956),  the 
Dirichlet  to  Neumann  (DtN)  boundary  condition  (Hunt,  1974),  infinite  elements  (Bettess,  1977), 
and  absorbing  boundary  layers  (Berenger,  1994)  are  all  techniques  that  seeks  to  eliminate  this 
reflection.  In  the  above  derivation,  a  simplification  of  the  Sommerfeld  radiation  condition  was 
used.  It  is  not  the  ideal  choice  for  a  radiation  condition  and  can  introduce  errors  into  the 
calculation.  However,  as  is  seen  in  a  later  section,  this  radiation  condition  is  necessary  to  solve 
the  problem  with  COMSOL  Multiphysics. 


3.3  Technique  for  Reducing  Numerical  Error  in  Two  and 
Three  Dimensions 

Increasing  the  number  of  elements  and  thus  nodes  in  the  solution  process  is  one  way  to 
reduce  the  error  in  a  finite  element  method  approximation.  Another  method  seeks  to  reduce  the 
random  numerical  error  after  the  solution  process  has  taken  place.  This  is  possible  in  some 
applications  of  acoustic  scattering  due  to  the  importance  of  point  measurements.  When 
determining  an  object's  target  strength,  it  is  not  necessary  to  solve  for  the  scattered  fields  at  all 
points.  Equations  (2.9)  and  (2.11)  show  how  only  a  single  point  measurement  of  scattered 
pressure  is  needed.  The  finite  element  method  approximates  the  solution  at  many  locations, 
namely  each  node  point.  If  the  approximation  is  examined  at  one  single  node  is  used,  there  will 
be  some  numerical  error  present.  However,  if  a  great  number  of  node  point  approximations  can 
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be  interpreted  to  determine  a  single  point's  value,  this  numerical  error  can  be  diminished. 
Techniques  are  now  presented  where  a  single  value  of  scattered  pressure  at  any  point  can  be 
calculated  from  the  knowledge  of  the  field  on  a  circle  or  a  sphere  surrounding  the  scatterer,  in 
two  and  three  dimensions  respectively.  The  location  of  the  calculated  pressure  is  arbitrary 
outside  of  the  immediate  vicinity  of  the  scatterer.  This  allows  the  scattered  pressure  to  be 
determined  at  points  well  beyond  the  finite  element  method's  numerical  domain.  In  this  way, 
scattered  pressure  can  be  determined  in  the  far-field  and  used  to  calculate  target  strength. 


3.3.1  Circular  Integral  Method 


First,  a  polar  coordinate  system  centered  at  the  scatterer  is  defined  with  radial  variable 
r  and  polar  variable  (p .  Start  with  the  Helmholtz  equation  for  the  scattered  field. 


V  p  ,  +  kzp  =0. 

r  scat  r  scat 


(3.23) 


The  Laplacian  operator  in  polar  coordinates  is  defined  as 


2  _  1  3  ,  3  .  1  3 

v  (r— )+— — 

r  or  or  r  o(p 


2  - 


(3.24) 


The  Helmholtz  equation  is  expanded  to  include  the  Laplacian, 


-T(r^fSL)+^^L+klp^=<>- 

r  or  dr  r  dcp 


(3.25) 


Assume  pscat  is  the  sum  of  a  unknown  number  of  separable  functions  of  r  and  (p , 


P  scat  ( r,g>)  =  'YjaT(r)PT((p) . 


(3.26) 


A  term  in  this  series  for  pscat  is  inserted  back  into  equation  (3.25),  and  separation  of  variables  is 
performed: 
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(3.27) 


1  d2j3r 


r  d  ,  dar,  ,  2  2 


2  2  c- 

(r-±)-k  r  =  ST 


PT  d(p~  az  dr  dr 

where  8n  is  a  constant.  Because  of  the  requirement  that  j3T(tp)  =  j8r{(p+  2 7t) ,  the  equation  in  (p 


can  be  solved  as 


P,=C,e 


in(p 


(3.28) 


with 


ST  =-n2, 


(3.29) 


where  n  is  any  integer.  The  expression  for  r  in  equation  (3.27)  becomes 

o  d2ar  dar 


+  r^L  +  aT(k2r2  —n2)  =0. 


(3.30) 


dr2  dr 


Performing  a  substitution  of  variables,  £  =  kr ,  this  becomes  the  well-known  equation 


e2  d  CC  yd(XT  .  e2  2\ 

^  37T- +  +  ar(¥ -n  )  =  0, 

d$  d% 


(3.31) 


whose  solutions  are  the  Bessel  functions  of  the  first  and  second  kind,  / „(£)  and  Y n ( g) 
(Abramowitz  and  Stegun,  1965).  Two  linear  independent  combinations  of  these  solutions  are 
the  Hankel  functions  of  the  first  and  second  kind,  7/ )''(£)  =  Jn(£)  +  ’ Yn ( %)  and 


a,  =  E,H<n(*r)  +  F,H“’(tr), 


(2), 


(3.32) 


where  E;  and  Fi  are  constants.  Because  (Xz  and  j3r  are  functions  of  integer  n  ,  the  summation 


in  equation  becomes 


(3.33) 


where  A,  and  B  are  constants. 
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The  n-dimensional  Sommerfeld  radiation  condition  (Sommerfeld,  1949)  is  given  as. 


lim  r  2  (--iku)  =  0. 

r^co  I  (lr 


( 3.34 ) 


For  the  two  dimensional  scattered  pressure,  the  condition  becomes 


lim  yfr 


ilcp  .  =  0. 


(3.35) 


Using  the  definition  of  the  derivatives  of  the  Flankel  functions  (Abramowitz  and  Stegun,  1965), 


H^\kr)'  =  ^-H^\kr)-HZikr) 


(3.36) 


H^\kr)'  =  ^-H^\kr)-H^(kr), 
kr 


the  term  in  equation  (3.35)  is 


.  -  An(—H"\kr)  -  HZ(kr)  -  iH?\kr)) 

_  ikPscat  )=Y^k^ 

+5„(f  H°\kr)  -  H^ikr)  - iH^\kr)) 
kr 


(3.37) 


Taking  the  limit  as  r  — >  °° , 


,, ,  2  i(Ar— — ) 

H„\kr )  — >  ,  —e  2  4 
V  7Tkr 


(3.38) 


2. ,  >  7l7t  71 . 
—i(kr - ) 

^r\kr)  — >  ,  ——e  24 


(3.39) 


li Bn(- 2  — e"  r_2_4  )  ** 

dr  V  # 


rw  . , ,  3/T . 

2/t  - -)s 


(3.40) 


To  satisfy  the  Sommerfeld  condition  in  equation  (3.35),  Bjn  =  0 .  The  expression  for  scattered 
pressure  becomes 
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(3.41) 


n=-oo 

Now  consider  a  circle  of  radius  centered  at  the  origin.  If  this  circle  completely  surrounds  the 
scatterer,  then  scattered  pressure  can  be  defined  on  it: 

=»  (3  42) 

n=-oo 

Multiplying  both  sides  by  e~im(p  and  integrating  over  (p  around  the  circle. 


l7t 


In 

J  pscat(R,(p)e-^Rd(p=  J  X  AnH"\kR)ein*e-im<pRcKp 


0  n= 


(3.43) 


In 


J  P„JR.<P)^"irdf>  =  2 zA.H“>(kR) 


(3.44) 


Solving  for  the  coefficient  An , 
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A 


l  l 

2 nhf(M)  0 


J  PscaAR,<P)e~,n(pd(p. 


(3.45) 


Therefore,  the  final  expression  for  the  scattered  pressure  at  point  (r,  <p)  depends  on  an  integral 
around  a  circle  of  radius  R  : 


In 


p~™ = htJiy me'"r  I  ■ 


(3.46) 


3.3.2  Spherical  Integral  Method 

In  three  dimensions,  consider  a  scatterer  located  at  the  origin  of  at  the  origin  of  a 
(r,  0 ,  (p )  spherical  coordinate  system.  Again  start  with  the  Helmholtz  equation,  (3.23).  The 
Laplacian  operator  in  spherical  coordinates  is  defined  as 
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(3.47) 


■nl  1  d  .  2  d  .  1  d  .  .  „  d  .  1  d 

v~  =  (r  4  +  ^TT77(sm  0—)  + 


rdr  dr'  r  sin#d#  dO  r2  smO  dcp1 


Assume  pscat  is  the  sum  of  a  number  of  separable  functions  of  r,  #,  and  (p : 


PscaM’0’<P)  =  5>r(r)A(0)n(P)- 


(3.48) 


A  term  in  this  series  is  inserted  back  into  the  Helmholtz  equation,  and  separation  of  variables  is 
performed. 


1  d2rr  sin2#  d  2dar.  sin#  d  „dB..  l2  2  •  2^, 

--4  = - —  (r  — 4 - ——(Sm0-ff)  -  42  sin2  0  =  et 

yT  d(p  aT  dr  dr  dO  dO 


(3.49) 


where  £t  is  a  constant.  Because  of  the  condition  yt((p)  =  yJ<p  +  2 71)  ,  the  equation  in  (p  can  be 


solved  as 


yr  =  cTel’\ 


(3.50) 


where  CT  is  a  constant  and  j  is  any  integer,  with 


£r=~J  ■ 

Inserting  (3.51)  into  (3.49),  separation  of  variables  is  performed  again, 

1  5  :(sin^)  -  - 4  -  - -4  ^L)  -*  V  =  8, . 


(3.51) 


(3.52) 


j8zsin0d0  d6  sin2#  dr  dr 


where  <?ris  a  constant.  Schwinger  et  al.  (1998)  give  that  the  equation  in  r  is  equal  to  the  product 


8  =  -n(n  + 1), 


(3.53) 


where  n  is  an  integer.  Substituting  x  =  cos  #  into  the  equation  for  #  gives  the  well  known 
equation 


(1-4 


9 


da: 


2  2.r^  +  A(4  +  l)-44=0, 


(3.54) 


l-.r 
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whose  solutions  are  the  Legendre  functions  of  the  first  and  second  kind  with  degree  n  and 
order  j ,  Pj(x)  and  Q,',(x)  (Abramowitz  and  Stegun,  1965).  Because  x  is  only  defined  on  the 
interval  [-1,1] ,  QJn  (x)  is  undefined  and  Pj(x)  =  PJ( cos  9)  is  the  only  solution.  Therefore, 

j3t  =  DrPnJ  (cos  9),  (3-55) 

where  Dt  is  a  constant.  The  equation  for  r  in  (3.52)  becomes 


— — —  (r2  — ‘-)  +  k2r 2  —n(n  + 1)  =0. 
CtL  dr  dr 


(3.56) 


A  substitution  g  =  kr  gives 


?^%-  +  2^  +  aT(?-n(n  +  l))  =  0, 


(3.57) 


whose  solutions  of  this  differential  equation  are  the  spherical  Bessel  functions  of  the  first  and 
second  kind,  j„(£)  and  (Abramowitz  and  Stegun,  1965).  Two  linear  independent 

combinations  of  these  are  the  spherical  Hankel  functions,  h(n1\^)  =  jn(£)  +  iy„(£)  and 


hin)(£)  = 


ar  =  Eftl'Hkr)  +  FThi2)(kr) , 


(3.58) 


with  constants  ET  and  Fz .  Because  (Xz ,  /3Z,  and  yz  are  functions  of  integers  j  and  n ,  the 
summation  for  scattered  pressure  in  equation  (3.48)  becomes 

°°  (3  59) 

P„„(r,e,<p)=  X  yi(At,h?(kr)  +  Bi„h?'(kr))p;( cos#)e'», 

"=-■> °  r=-“ 

where  A  jn  and  Bjn  are  constants.  Recall  the  three  dimensional  Sommerfeld  radiation  condition 
in  equation  (3.6).  Using  the  definition  of  the  derivatives  of  spherical  Hankel  functions. 
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(3.60) 


K\kr)'  =  ^\kr)-hZ{kr) 
kr 

K2\kry  =  ^h?\kr)-^(kr), 

kr 

the  term  in  the  limit  of  the  Sommerfeld  condition  becomes 

r(<!EsL-ikPmi)  = 

dr 


(3.61) 


"  h"\kr)  -  h™(kr)  -  ih™(kr )] 

kr 

+Bjn[^hl2\kr)  -  lC(kr)—ihl2\kr)] 
V  kr 


PJ  (cos  9)elJf . 


Taking  the  limit  r  — >  °° , 


l  i(kr™-JL) 

li^ikr  )  — > — e  2  2 

kr 

(3.62) 

h(n2\kr)  — >  — e  2  2 

kr 

(3.63) 

oo  oo  (  ...  n7t  ^ 

P^(cosd)e"f 

LU  j=-oon=-oo\  y 

(3.64) 

To  satisfy  the  Sommerfeld  condition,  B  jn  =  0 .  Equation  (3.59)  becomes 

PscaM’di<P)  =  X  Z^CWACOS^. 

«=-°°  7=-°° 

(3.65) 

Legendre  functions  follow  two  important  recurrence  relations  (Abramowitz 

and  Stegun,  1965): 

^+1(z)  =  (z2  - 1)-1/2  [(v  -  M)zPy(z)  -  (v  +  ju)PUz )] 

(3.66) 

ii 

(3.67) 
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These  allow  any  Legendre  function  Pr'  to  be  written  as  a  combination  of  others  with  both  n 
positive  and  j  e  [— n,  —  n  + 1 . . .  n  —  1,  n] .  This  simplifies  the  summation  limits  in  the  expression  for 
scattered  pressure: 


=  zt  AiX"(kr)Pj( cos»)e‘Ir 

n= 0  j=—n 


(3.68) 


Now  consider  a  sphere  of  radius  R  centered  at  the  origin  with  coordinates  (R.  0o,<po)  .  If  this 


completely  surrounds  the  scatterer,  the  scattered  pressure  can  be  defined  on  it: 

=  tt  <c°s • 

n= 0  j=—n 

Multiply  both  sides  by  Pj( cos  d0)e~ij<Po ,  where  n  and  j  are  each  some  integer, 

Psca,(R’  0(y<P,)Pj(COS  0{))e~ii%  = 


(3.69) 


(3.70) 


X  X  AjnKl)(kR)Pj  (cos  &0  )e>i%Pj  (cos  0O  )e~'J 


Wo 


n= 0  j=—n 


This  expression  is  then  integrated  over  the  surface  of  the  sphere, 
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{  {  Psca, (R,  *„,#>)# (cos d0)e-™R2  sin  d0d%ddQ 


0  0 


(3.71) 


n  In 


JILL  Ajnh(n[ ' ( kR ) PJ (cos  00 )PJ (cos  0O )e'I  J  J,'P"R2  sin  0()d %d 6{) 


0  0  n=0  j=~n 


The  integral  over  (p()  on  the  right  is  zero  when  j  ^  j  and  2n  when  j  =  j .  Therefore  the 
expression  simplifies  to 


n  In 

1 1  p,rJRA,n)Ph^e„)e-ilK 


0  0 


sin  00d  <p()d  60 


(3.72) 
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lK\  X  Ajnh«]  (kR)P,!  (COS  °0  )Pi  (C0S  P0  )  sin  60d P0  • 

0  n=0 


A  substitution  of  variables  cos  0O  =  a:  into  the  right  side  of  the  expression  gives 


n  2  n 


J  J  (cos^V^ sin 60d<p0dd0  =  2n\ £ A,hl1\kR)PJ(x)PJ(x)dx . 

_1  n=0 


(3.73) 


0  0 


Legendre  functions  of  the  same  degree  and  different  order  are  orthogonal  on  the  interval  [—1,1] 
(Abramowitz  and  Stegun,  1965): 

(3.74) 


J  p:(x)p;‘(x)dx  =  o, 


-1 


when  b  =£  c .  Therefore,  the  expression  becomes 
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J  J  Psca, ( R>  > % ) pi! (cos d<M  m  Sin  d0d<p0dd0  =  A,M£\kR)\ Pj(x)2dx  . 


(3.75) 


0  0 


Now  A.  can  be  solved  for. 


jn 
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J  J  P^iR, ^<P0)p,!(cos d0)e~^  sin d,d%dd. 


(3.76) 


A..  = 


_  00 


jn 


i 

2  7Th™(kR)\ Pj(x)2dx 


-1 


Two  cases  must  be  considered:  j  >  0  and  j  <  0  .  When  j  >  0 ,  there  is  an  expression  for  the 


integral  in  the  denominator  of  (3.76)  (Abramowitz  and  Stegun,  1965): 

(n  +  j) ! 


\PJ(x)2dx  = 


(3.77) 


(n  +  1/  2){n  -  j) ! 


Inserting  this  back  into  the  expressions  for  A.  and  scattered  pressure  in  (3.68), 
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J  J  Pscai  ( R ,  00 -  <Po)Pj  (cos  0O )e~m  sin  0od %d 0O  (n  +  l/2)(n-  j ) ! 


(3.78) 


4  =  Q  0 

jn 


2  jvh™(kR)(n  +  j)\ 


P,Car(r^^)  =  — ZZ 

n=0  j=i 


(/?  +  !/  2)(n  -  j )  cos 

(n  +  jy.h^m 


(3.79) 


<|  J  /L™, ( #o •  % ) ^'(cos 0O )U'Wl  sin  B{]d <pt)d 0O . 


The  second  case  is  when  j  <  0 .  Legendre  functions  of  negative  order  can  be  written  as 


pj(x)  =  (»  +  j  +  Wpi(x) 

(n-i+l)\ 


(3.80) 


Equation  (3.76)  can  be  rewritten  as 


J  J  VscaAR^0 o,<po)Pj(cos0Q)e~ij%  sin 0od<pod6o 


A  -  0  0 

jn 


2 Jcti£\kR)  f  \ {n  +  j  +  1)l  P;j(x)  1  dx 


(3.81) 


Now  the  integral  in  the  denominator  can  be  evaluated: 


J  J  PsaJP- 0o><Po )Pj(cos 0o)e-m  sin 60d(p0d 60 


A  -  0  0 

jn 


2  Kh^ikR) 


(n  +  j  +  l)\\  (n  -  j) ! 


(n-j  +  l)\  (n  + 1  /  2)(n  +  j) ! 


(3.82) 


Plugging  this  back  into  the  expression  for  scattered  pressure  (3.68) , 


,.0n,_  1  (n  +  l/2)(n  +  j)\hV\kr)Pj(cos0)eijr 

P scaAr >rV  /  \  2 

(„-;)!  <^ W) 

{(n-j  +  l)\) 


(3.83) 


<J  J  P-atiR,  Polices  0o)e~^  sin  0Qd(podOo 


Using  the  expression  in  (3.80),  this  becomes 
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(3.84) 


PscaM^^)  =  —YjYj 

n—O  j=—n 


1  ^  ^  in  + 1  / 2 ){n  +  j)\h™(kr)P;J(c os 0)^ 


(n-j)\h?(kR) 


7t  27t 


XJ  J  PsrJR.e{],%)P„  J(cos0ft)e^m'  sin 60d<p0dd0 . 


0  0 


Comparing  (3.79)  when  j  >  0  to  (3.84)  when  j  <  0  ,  it  is  clear  the  expression  for  scattered 
pressure  may  be  written  in  one  form  for  all  j : 


<  0  ,  1  v  V(n  +  1/2)(H  j\y-^\kr)P}!\cosdyj<p 

P scat\r >rv  2  /  -  /  - 


n=0  j=—n 


(n  +  \j\VXl\kR) 


(3.85) 
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XJ  J  PscaSR’  (cos£0)c~'mj  sin  6()d(p0dd{) 


0  0 


This  equation  can  be  written  so  that  the  surface  integral  is  more  general: 

°°  n  U^(br\ 

p_('.W  =  EZ 

,1=0  j=-« \KK)  r 


(3.86) 


and 


<$>jn(R,9,(p) 


|(2»  +  !)(»-  \p 

]j  47rR2(n  +  |j|)!  " 


(cos  0)e‘1,p , 


where  T  is  the  surface  of  the  sphere  and  *  denotes  the  complex  conjugate. 

Equations  (3.46)  and  (3.86)  give  the  expressions  for  scattered  pressure  at  any  point  in 
two  and  three  dimensions,  based  on  the  knowledge  of  the  field  on  a  circle  and  a  sphere 
surrounding  the  scatterer  respectively.  These  expressions  are  very  helpful  in  conjunction  with 
the  finite  element  method,  which  solves  for  the  scattered  pressure  over  a  finite  domain  around 
the  object.  The  scattered  pressure  can  be  ascertained  on  a  circle  or  sphere  from  the  finite  element 
method,  and  then  used  to  find  its  value  at  any  point  by  the  above  techniques.  In  addition,  they 
prove  to  help  remove  significant  amounts  of  numerical  error. 
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3.4  Implementation  of  COMSOL  Multiphysics 


COMSOL  Multiphysics  (ver.  3.2a)  was  used  in  this  work  for  the  finite  element  method 
simulations. 

3.4.1  COMSOL  Methodology 

3.4.1.1  Geometry  and  Elements 

The  process  of  creating  scatterer  geometry  within  COMSOL  was  discussed  when  solving 
the  Kirchhoff  integral  numerically.  In  finite  element  method  scattering  problems,  the 
surrounding  medium  subdomain  must  be  created  as  well.  To  simplify  the  implementation  of 
the  radiation  condition  on  the  boundary,  this  exterior  region  is  modeled  as  a  sphere  in  three 
dimensions  and  a  circle  in  two.  The  radius  of  this  sphere  (or  circle)  is  labeled  as  A .  To  help 
visualize,  figure  28  shows  two  sample  finite  element  domains  in  COMSOL.  The  sphere  and 
circle  represent  the  surrounding  medium  subdomains,  while  the  box  and  square  represent  the 
scatterer  subdomains. 
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Figure  28  -  Sample  three  dimensional  and  two  dimensional  finite  element  domains  in  COMSOL. 

After  the  geometries  of  these  two  subdomains  —  scatterer  and  surrounding  medium— are 
defined,  the  elements  themselves  can  be  created.  Identical  to  the  process  in  the  Kirchhoff 
chapter,  this  consists  of  discretizing  the  domain  into  quadrilaterals  (or  triangles  in  two 
dimensions).  In  this  application,  the  size  of  the  elements  is  again  important  to  the  accuracy  of 
the  solution.  Elements  may  be  defined  by  a  maximum  allowable  size,  /imax  .  The  effects  of 
element  size  will  be  explored  shortly. 


3.4.1.2  Defining  the  Partial  Differential  Equations 

The  next  step  is  to  create  the  partial  differential  equations  that  are  used  in  the  numerical 
calculation.  This  is  done  by  defining  the  coefficients  of  the  generic  partial  differential  equations 
used  in  COMSOL.  Every  subdomain  and  boundary  between  subdomains  requires  an  equation. 
In  the  case  of  acoustic  scattering,  the  subdomain  equations  describe  the  traveling  waves  in 
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terms  of  total  pressure,  the  sum  of  incident  and  scattered  fields.  Boundary  equations  are  used  to 
represent  both  the  boundary  condition  between  the  scatterer  and  surrounding  medium  and  the 
radiation  condition  at  the  edge  of  the  numerical  domain. 

The  general  equation  for  subdomains  in  COMSOL  is 

V  •  (— cVu  -  OCu  +  y)  +  an  +  j3  ■  Vu  =  / ,  (3.87) 


where  the  coefficients  must  be  defined.  Let  u  represent  the  total  pressure: 


U  =  Pro,  =  Pine  +  Psca,  • 


(3.88) 


From  COMSOL' s  solution  for  total  pressure,  scattered  pressure  is  easily  solved  for  by 
subtracting  out  the  known  incident  field.  To  describe  the  pressure  field,  equation  (3.87)  should 
model  the  Helmholtz  equation  (3.1).  This  is  easily  done  by  configuring  the  coefficients  as 
follows: 


c  =  1 
a  =  0 
y  =  0 

a  =  - k 2 

fi  =  0 
f  =  0 


(3.89) 


Terms  like  k  can  either  have  their  values  entered  explicitly  into  the  equations,  or  they  can  be 
defined  separately  to  allow  for  symbolic  equations. 

The  first  boundary  condition  that  has  to  be  described  is  between  the  scatterer  and  the 
surrounding  medium.  This  can  be  modeled  as  a  variety  of  different  boundaries:  rigid/ fixed, 
soft,  fluid,  and  potentially  elastic  (although  this  was  not  tested).  COMSOL  provides  two  generic 
equations  that  can  be  used  to  model  a  Neumann  or  Dirichlet  boundary  condition, 

n-(cVu  +  au-  y)  +  qu  =  g  -hT n  (3.90) 
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hu  =  r , 


where  n  is  the  outward  normal  vector  from  the  surface  and  /U  is  a  Lagrange  multiplier.  Several 
of  the  coefficients  here  (c,CC,y)  are  already  defined  in  the  subdomain  description  in  (3.89). 
Substituting  these  values  and  rewriting, 

n-Vu  + qu  =  g  -  hJ  jU  (3.91) 

hu  =  r . 

By  manipulating  the  coefficients  of  this  equation,  all  of  the  necessary  boundary  conditions  can 
be  replicated.  For  example,  in  a  rigid/ fixed  boundary,  the  normal  gradient  of  pressure  is  zero 
on  the  boundary.  To  model  this,  a  Neumann  boundary  condition  is  used  and  all  remaining 
coefficients  are  set  to  zero: 

n-Vu  =  0  (3.92) 

0  =  0. 

For  a  soft  boundary  condition,  the  pressure  goes  to  zero  at  the  interface.  To  model  this,  h  is  set 
at  a  nonzero  value  in  order  to  include  the  second  equation.  The  simplest  method  is  to  set  it 
equal  to  one.  Then  because  of  the  Lagrange  multiplier,  the  first  equation  imposes  no  restriction 
on  u  and  can  thereby  be  ignored.  The  second  equation  simply  ends  up  being 

u  =  r  (3.93) 

and  r  can  be  set  to  any  value;  this  is  zero  in  the  case  of  the  soft  boundary  condition.  Similar  steps 
can  be  taken  to  create  a  fluid  boundary,  requiring  continuity  of  both  pressure  and  pressure 
gradient  at  the  interface.  In  this  work,  the  rigid/fixed  condition  was  always  used. 

The  need  to  model  free  space  at  the  boundary  of  the  finite  surrounding  medium 
subdomain  in  acoustic  scattering  problems  has  been  discussed.  One  way  to  do  this  is  to  use  a 
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simplified  version  of  the  Sommerfeld  radiation  condition,  introduced  in  a  previous  section  as 
equation  (3.7), 


dPscat 

dr 


=  ikp 


scat  * 


(3.94) 


This  is  an  equation  for  the  scattered  field,  not  the  total  field  for  which  COMSOL  solves. 
Therefore,  the  scattered  field  must  be  written  in  terms  of  the  total  field  and  the  known  incident 
field. 


Psca,  =  U  ~  P, 


(3.95) 


The  simplified  Sommerfeld  condition  becomes 

77  ■  V(w  -  pmc)  =  ik(u  -  pmc)  . 


(3.96) 


The  incident  field  must  now  be  defined.  For  an  incident  plane  wave  of  amplitude  pQ  and  wave 


(3.97) 

(3.98) 


Looking  again  at  COMSOL' s  equations  in  (3.90),  this  condition  can  be  represented  by 

configuring  the  coefficients  as  follows: 

q  =  -ik  (3.99) 

g  =  i  (n-k)-k  p0e‘^'r) 
h  =  0 

r  -  0 
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3.4.1.3  The  Solution  Phase 


The  large  matrix  equations  are  automatically  compiled  by  COMSOL  and  run  through  one 
of  its  linear  solvers.  For  two-dimensional  problems,  the  UMFPACK  solver  was  used.  This  solves 
the  matrix  equations  directly  using  LU  factorization.  For  the  three-dimensional  problems 
requiring  much  more  memory,  the  GMRES  iterative  solver  was  used  instead.  While  a  slower 
solver  than  UMFPACK,  GMRES  allowed  for  calculations  with  many  more  elements. 


3.4.2  Accuracy  of  the  Solution 


The  accuracy  of  the  finite  element  method  approximation  can  be  judged  several  ways.  It 
can  simply  be  interpreted  at  a  single  node  point  as  the  difference  between  the  exact  solution  and 
the  computed  value. 


e(x)  =  u(x)-uFEM(x). 


(3.100) 


There  are  a  variety  of  error  norms  that  give  a  more  global  understanding  of  a  solution's 
accuracy.  Two  examples  are  the  energy  and  maximum  norms  (Becker  et  ah,  1981), 

1 1/2  (3.101) 


||  e\\E=\  | [(e(x)')2  +  e(x)2~\d£l  l  , 


I  e  11^=  max  I  e(x)  I. 


(3.102) 


Whenever  using  the  finite  element  method  (or  any  numerical  technique),  error  in  the 
solution  is  expected  as  a  result  of  simply  discretizing  the  domain.  Other  error  effects  are  specific 
to  the  problem;  in  this  case,  it  is  the  use  of  the  simplified  Sommerfeld  radiation  condition.  In 
order  to  help  understand  and  attempt  to  separate  the  effects  from  the  discretization  and  the 
radiation  condition,  two  unitless  parameters  were  created:  EPW  and  SMD . 
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Given  a  mesh  of  elements  with  maximum  element  size  hm.a  for  a  problem  with 
wavenumber  k , 


EPW  = 


X 


h 


2  n 

h,„k 


(3.103) 


In  this  way,  there  are  at  least  EPW  elements  per  wavelength.  By  increasing  EPW  ,  the 
elements  will  get  smaller  relative  to  wavelength.  Therefore  this  parameter  is  a  measure  of  the 
discretization,  generally  with  larger  values  resulting  in  more  accurate  solutions. 

For  a  spherical  surrounding  medium  domain  of  radius  A  centered  at  the  origin,  with  a 
scatterer  whose  maximum  extension  from  the  origin  is  distance  /ax ,  and  with  wavenumber  k , 
define 


SMD  =  k{A-rmJ. 


(3.104) 


This  parameter  is  the  minimum  distance  a  scattered  wave  must  travel  before  it  reaches  the  outer 
boundary  of  the  surrounding  medium  in  terms  of  wavenumber.  Because  the  simplified 


Sommerfeld  radiation  condition  has  an  error  of  order  o 


rn 

U; 


,  increasing  SMD  should 


increasing  the  accuracy  of  the  solution. 


3.4.2.1  Two  Dimensional  Test:  Rigid/Fixed  Infinite  Cylinder  (Max  Error  Norm) 

The  finite  element  method  with  the  simplified  Sommerfeld  condition  can  be  tested  for 
several  problems  with  known  solutions.  A  simple  problem  that  can  be  tested  is  scattering  from 
a  rigid/ fixed  infinite  cylinder,  which  has  a  known  scattering  solution  that  does  not  vary  along 
its  length  (Rayleigh,  1945).  For  an  incident  plane  wave  pmc  =  p0e,kx , 


84 


(3.105) 


and 


Psca,  =  Po2X  cos  (m<p)H^(kr) 

m—0 

—£mi"‘ J  m(ka)' 


(3.106) 


where  £m  is  equal  to  1  for  m  =  0  and  equal  to  2  for  all  other  vales  of  m  .  Because  of  the 

symmetry  of  the  infinite  cylinder,  this  problem  may  be  thought  of  as  the  scattering  from  a 
rigid/ fixed  circle  in  two  dimensions. 

The  two  dimensional  problem  can  be  run  in  COMSOL,  and  the  approximations  at  each 
node  can  be  compared  to  the  exact  solutions  at  those  points.  Assume  an  infinite  cylinder  lies 
along  the  z  axis  with  radius  a  .  The  cylinder  experiences  an  incident  field  of  e,kx ,  representing  a 
plane  wave  traveling  in  the  positive  x  direction  with  amplitude  p0  =  1  and  wavenumber  k  .  To 

model  the  problem  in  COMSOL,  a  circle  of  radius  a  =  1  centered  at  the  origin  represents  the 
cylinder.  Because  of  symmetry, 

SMD  =  k(A-a) .  (3.107) 

This  parameter  was  varied  over  several  values  and  was  used  to  determine  A ,  the  radius  of  the 
outer  domain.  The  parameter  EPW  was  varied  from  two  to  eight,  representing  the  minimum 
number  of  elements  per  wavelength.  The  solution  was  generated  in  COMSOL  and  the  scattered 
field  was  determined  at  each  node  point. 
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Surface:  p_tot 


Surface:  p_scat 


Figure  29  -  Total  and  scattered  pressure  fields  from  a  rigid/fixed  infinite  cylinder,  ka  =  10  in 
COMSOL.  The  incident  wave  arrives  from  the  left  side  of  the  domain. 

The  maximum  error  norm  was  determined  by  comparing  the  solutions  at  every  node  and 
finding  the  maximum  deviation  from  the  exact  solution.  This  was  done  at  two  separate 
frequencies,  determined  by  ka,  the  product  of  the  wavenumber  and  the  cylinder's  radius. 
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Max  error  norm  for  a  rigid/fixed  infinite  cylinder,  ka  =  1 
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Figure  30  -  Max  error  norm  versus  EPW  for  COMSOL  approximation  of  scattering  from  a  rigid/fixed 
infinite  cylinder  at  ka  =  1. 
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Figure  31  -  Max  error  norm  versus  EPW  for  COMSOL  approximation  of  scattering  from  a  rigid/fixed 
infinite  cylinder  at  ka  =  5. 
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As  expected,  decreasing  element  size  also  decreases  the  error  norm  of  the  solution.  The 
error  slows  its  reduction  and  in  some  cases  stops  it  as  when  the  elements  reach  a  certain  size. 
The  error  that  remains  for  the  curves  that  have  leveled  out  is  a  result  of  the  radiation  condition 
used  and  cannot  be  fixed  by  discretizing  the  domain  further.  As  expected,  within  this  range  of 
small  element  size,  domain  size  is  a  key  factor  in  the  amount  of  error  in  the  solution.  Increasing 
the  domain  size  helps  to  reduce  error  when  further  discretization  has  no  effect. 


3A.2.2  Two  Dimensional  Test:  Rigid/Fixed  Infinite  Cylinder  (Infinite  Form  Function) 


Because  target  strength  is  undefined  for  infinite  objects,  the  infinite  form  function  /” 
can  be  used  for  comparison  instead.  Similar  to  the  scattering  amplitude  of  a  finite  object,  the 
infinite  form  function  of  a  cylinder  of  radius  a  is  defined  only  in  the  far-field  kr  »  1  (Stanton, 
1992): 


^  ikr  r  °° 

P  scat  *  f  ' 


(3.108) 


For  a  cylinder  aligned  with  the  z  axis,  the  infinite  form  function  for  a  rigid/fixed  cylinder  from 


an  incident  plane  wave  pmc  =  p0e,kx  is  defined  as 


r  = 


-=em'AYJbm  cos (m<p) 
v  71  ka  m=o 


(3.109) 


and 


b,„  = 


-£mJm(ka) 
H(:\ka)'  ' 


(3.110) 


In  the  backscattering  direction,  (p  =  7l ,  giving 


fbs  = 


-in!  A 


V  71  ka 


IX(-D” 


(3.111) 


m—0 


Again,  the  only  interest  here  is  in  this  backscattering  direction,  so  the  subscript  will  be  dropped. 


88 


The  value  of  scattered  pressure  used  to  find  the  infinite  form  function  in  equation  (3.108) 
can  be  determined  two  ways  from  COMSOL's  solution.  It  can  either  be  taken  as  a  simple  point 
measurement  from  the  solution  itself,  or  it  can  be  calculated  using  the  circular  integral  method 
in  equation  (3.46).  This  latter  method  was  used  by  taking  the  solution  at  a  series  of  points  on  a 
circle  surrounding  the  scatterer. 

If  pscat  is  taken  as  a  point  measurement,  in  order  to  ensure  that  it  is  in  the  far-field,  it  is 

determined  at  distance  d  from  the  scatterer  where  kd  »  1 .  In  the  example  below,  kd  =  20  . 

The  size  of  the  domain  was  based  on  SMD  =  50 ,  and  element  size  was  based  on  EPW  =  6  . 


Figure  32  -  Absolute  value  of  infinite  form  function  of  a  rigid/fixed  infinite  cylinder,  using  the  modal 
series  solution  and  COMSOL  (point  method).  EPW  =  6  and  SMD  =  50. 

There  are  errors  throughout  most  of  the  ka  region,  with  them  growing  with  ka. 
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For  the  circular  integral  method,  the  same  values  for  SMD  and  EPW  were  used.  The 
integral  was  performed  on  a  circle  with  radius  1.1  times  the  cylinder  radius,  surrounding  the 
scatterer. 
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Figure  33  -  Absolute  value  of  infinite  form  function  of  a  rigid/fixed  infinite  cylinder,  using  the  modal 
series  solution  and  COMSOL  (circular  integral  method).  EPW  =  6  and  SMD  =  50. 


The  circular  integral  method  succeeded  in  removing  much  of  the  error  from  the  point 
measurement  method.  There  is  still  some  small  error  visible,  caused  by  the  discretization  and 
radiation  condition. 

To  determine  whether  or  not  the  solution  has  converged  with  the  given  values  of  SMD 
and  EPW  ,  convergence  tests  were  run.  First,  setting  SMD  =  50 ,  and  varying  EPW  ,  the 
infinite  form  function  was  found  for  different  values  of  ka  using  the  circular  integral  method. 
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Infinite  form  function  error  for  a  rigid/fixed  infinite  cylinder,  SMD  =  50 
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Figure  34  -  Error  in  the  infinite  form  function  approximation  of  a  rigid/fixed  infinite  cylinder  at  SMD 
=  50. 


The  solution  has  definitely  converged  to  a  very  small  error  at  EPW  >  6  for  both  cases.  Next 
setting  EPW  =  6 ,  the  same  process  was  run  while  varying  SMD . 
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Infinite  form  function  error  for  a  rigid/fixed  infinite  cylinder,  EPW  =  6 
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Figure  35  -  Error  in  the  infinite  form  function  approximation  of  a  rigid/fixed  infinite  cylinder  at  EPW 

=  6. 


It  first  should  be  noted  that  the  error  scale  used  in  this  figure  is  much  smaller  than  that  of  the 
previous  figure.  Therefore,  although  the  approximations  seem  to  oscillate  at  higher  values  of 
SMD ,  the  solution  seems  to  have  converged  into  a  very  small  error  region.  From  this  figure  it  is 
also  seen  that  at  a  certain  point,  increasing  domain  size  does  not  seem  to  reduce  the  error  effects 
of  the  radiation  condition. 


3A.2.3  Three  Dimensional  Test:  Rigid/Fixed  Sphere 

A  three  dimensional  scattering  problem  with  known  solution  is  that  of  a  rigid/ fixed 
sphere  (equation  (2.19)).  The  added  dimension  greatly  increases  the  number  of  finite  elements 


92 


for  a  given  domain  size.  Often  the  domain  of  the  problem  cannot  be  extended  into  the  far-field 
due  to  computational  constraints.  Therefore,  the  spherical  integral  method  proves  very  useful  in 
allowing  near-field  data  to  be  used  to  calculate  a  value  in  the  far-field  so  that  target  strength 
may  be  determined.  Available  computer  resources  limited  three  dimensional  COMSOL 
simulations  to  the  use  oiEPW  =  6  and  SMD  =  20,  with  each  simulation  requiring  several 
hours.  Using  these  parameters,  the  target  strength  of  a  rigid/fixed  sphere,  radius  5  cm  was 
calculated  and  compared  with  the  exact  modal  series  solution  (section  2.3.2). 
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Figure  36  -  Target  strength  of  a  rigid/fixed  sphere,  radius  5  cm,  using  the  modal  series  solution  and 
COMSOL  (spherical  integral  method).  EPW  =  6  and  SMD  =  20. 

At  values  of  ka  less  than  five,  the  finite  element  results  are  very  good.  They  are  able  to  capture 
the  proper  peak  null  structure  of  the  exact  solution  much  better  than  the  Kirchhoff  method.  This 
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is  because  of  the  finite  element  method  does  not  ignore  the  effects  of  diffraction.  However  at 
higher  values  of  ka,  the  results  begin  to  diverge  from  the  exact  solution,  causing  errors  on  the 
order  of  2-3  dB.  This  is  an  effect  of  some  error  in  the  finite  element  method,  most  likely  caused 
by  the  radiation  condition  used. 


3.4.3  Methodology  for  Sand  Dollar  Predictions 

The  process  of  predicting  scattering  from  any  three  dimensional  object  is  similar  to  that  of 
the  sphere  in  the  previous  section.  The  only  difference  lies  in  the  creation  of  the  scatterer 
geometry,  which  has  been  discussed.  For  the  predictions  of  the  scattering  from  the  sand  dollar, 
a  disk  was  used  to  simplify  this  geometry.  As  was  seen  with  the  rigid/ fixed  sphere,  errors  on 
the  order  of  2-3  dB  were  expected  from  the  use  of  COMSOL,  likely  stemming  from  the  choice  of 
radiation  condition.  However,  the  experimental  sand  dollar  scattering  data  were  such  that  these 
errors  were  small  enough  for  the  finite  element  method  to  be  useful.  The  results  are  presented 
and  compared  with  experimental  data  in  Chapter  4. 


3.5  Chapter  Summary 

The  finite  element  method  has  been  explored  in  the  application  of  acoustic  scattering 
problems.  Techniques  for  reducing  numerical  error  after  the  solution  process  in  two  and  three 
dimensions  were  presented.  These  techniques,  referred  to  here  as  the  circular  and  spherical 
integral  methods,  also  allow  scattered  pressure  values  to  be  calculated  at  points  outside  of  the 
original  domain.  COMSOL  Multiphysics  was  explored  in  the  simple  tests  of  scattering  from  a 
rigid/fixed  infinite  cylinder  and  a  sphere.  The  finite  element  method  is  used  to  predict  the 
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scattering  from  rigid/ fixed  disks  in  Chapter  4.  These  results  are  compared  to  the  experimental 
data  of  scattering  from  a  sand  dollar. 
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Chapter  4  Model  Results  and 
Comparison  to  Laboratory 
Experimental  Data 

4.1  Introduction 

In  this  chapter,  the  sand  dollar  scattering  models  are  explored,  and  predictions  are 
compared  to  laboratory  experimental  results.  The  models  consist  of  the  analytic  and  numerical 
solutions  of  the  Kirchhoff  integral  and  solutions  from  the  finite  element  method.  Before 
introducing  the  experimental  data,  the  model  predictions  are  compared  with  one  another.  All 
model  predictions  are  based  on  rigid/  fixed  boundary  conditions.  The  analytic  solution  of  the 
Kirchhoff  integral  for  a  disk,  chosen  as  a  simple  approximation  to  the  geometry  of  the  flat  side 
of  the  sand  dollar,  is  compared  with  the  predicted  results  from  the  finite  element  method  for  a 
disk  over  both  ranges  of  frequency  and  orientation.  Then,  the  model  predictions  for  the 
rigid/ fixed  disk  and  spherical  cap  are  compared  to  the  experimental  results  from  Stanton  and 
Chu  (2004)  of  the  scattering  from  an  aluminum  disk  and  a  sand  dollar.  These  comparisons  are 
conducted  over  a  range  of  angles  of  orientation,  -20  to  80  degrees  from  broadside,  at  a  single 
frequency,  70  kHz.  This  range  of  angles  was  chosen  due  to  symmetry  and  difficulties  in  the 
finite  element  method  near  90  degrees.  In  addition  to  the  analytic  solutions  of  the  Kirchhoff 
integral  and  the  finite  element  method  for  a  rigid/  fixed  disk,  a  numerical  solution  of  the 
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Kirchhoff  integral  is  also  compared  with  the  experimental  results.  This  numerical  solution  is 
found  for  a  rigid/fixed  three  dimensional  model  of  the  sand  dollar  from  the  experiments, 
assembled  from  CT  scans.  Finally,  there  is  significant  evidence  that  rigid/fixed  boundary 
conditions  are  not  appropriate  for  an  aluminum  disk,  and  a  heuristic  approach  for  accounting 
for  the  penetrable  boundary  condition  is  presented.  This  correction  is  also  applied  to  the  model 
predictions  of  scattering  from  the  sand  dollar,  as  the  boundary  conditions  appropriate  for  a 
sand  dollar  are  not  accurately  known. 


4.2  Comparison  of  Predicted  Scattering  Based  on  the 
Kirchhoff  and  Finite  Element  Methods  for  a 
Rigid/ Fixed  Disk 

In  this  section,  the  predictions  from  the  analytic  Kirchhoff  method  and  finite  element 
method  are  compared.  The  two  models  were  used  to  calculate  the  target  strength  of  a 
rigid/ fixed  disk  over  ranges  of  both  frequency  and  angle  of  orientation.  The  analytic  solution  of 
the  Kirchhoff  method  was  derived  for  a  rigid/ fixed  disk  in  section  2.5.2.  For  the  finite  element 
method  solution,  the  program  COMSOL  was  used.  The  disk  was  given  a  radius  of  3.625  cm  and 
a  thickness  of  5.5  mm.  This  size  was  chosen  based  on  the  actual  dimensions  of  the  sand  dollar 
from  Stanton  and  Chu  (2004),  which  will  be  seen  in  the  following  section. 

The  frequency  dependent  calculations  were  performed  with  the  incident  acoustic  wave  at 
broadside  to  the  disk,  allowing  the  finite  element  simulations  to  be  run  in  two  dimensions, 
resulting  in  a  much  faster  and  more  efficient  calculation  with  EPW  =  6  and  SMD  =  50 . 
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Afterwards,  the  spherical  integral  method  described  in  Chapter  3  was  used  to  calculate  target 
strength.  The  Kirchhoff  and  finite  element  methods  were  carried  out  over  a  range  of  ka,  a 
dimensionless  quantity  given  by  the  product  of  the  wavenumber  and  the  disk's  radius,  as 
illustrated  in  figure  37. 
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Figure  37  -  Target  strength  at  broadside  incidence  of  a  rigid/fixed  disk  of  radius  3.625  cm  and  thickness 
5.5  mm,  based  on  (1)  the  analytic  solution  of  the  Kirchhoff  integral  (dashed  line)  and  (2)  the  finite 
element  method  using  COMSOL  with  axial  symmetry  (line  with  x's).  EPW  =  6  and  SMD  =  50. 


There  is  very  good  agreement  over  most  of  the  range  of  ka  for  the  two  models.  The  undulations 
in  the  finite  element  method  predictions,  particularly  apparent  at  low  ka,  most  likely  result 
from  diffraction  effects.  These  are  not  accounted  for  in  the  Kirchhoff  method,  as  was  also  seen 
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when  the  results  were  compared  with  the  modal  series  based  solutions  for  the  sphere  and  finite 
cylinder  in  Chapter  2. 

The  target  strength  of  the  disk  was  also  calculated  at  a  single  frequency,  70  kHz,  and 
over  a  range  of  angles  of  orientation,  from  zero  (broadside)  to  eighty  degrees,  seen  in  figure  38. 
The  disk  was  rotated  in  increments  of  one  degree  around  an  axis  perpendicular  to  the  incident 
wave  vector.  Because  of  symmetry,  it  only  had  to  be  rotated  in  one  direction.  The  frequency  of 
70  kHz  was  chosen  to  match  that  of  the  experimental  data  in  the  following  section  and 
corresponds  to  ka  =  10.6. 


Disk  (Kirchhoff  and  finite  element  methods) 


Figure  38  -  Target  strength  at  70  kHz  of  a  rigid/fixed  disk  of  radius  3.625  cm  and  thickness  5.5  mm 
based  on  (1)  the  analytic  solution  of  the  Kirchhoff  integral  (dashed  line)  and  (2)  the  finite  element 
method  using  COMSOL  (line  with  x's).  EPW  =  6  and  SMD  =  20. 


99 


The  predictions  are  in  relatively  good  agreement  at  broadside  (angle  of  orientation  zero), 
expected  because  of  the  good  agreement  at  ka  =  10.6  in  the  frequency  dependent  comparison, 
though  the  finite  element  method  predicted  a  slightly  larger  (1.1  dB)  target  strength  at 
broadside.  The  side  lobe  structures  also  match  up  relatively  well,  particularly  in  the  location  of 
the  nulls,  although  there  are  height  differences  of  up  to  5  dB  in  target  strength  on  some  of  the 
peaks.  The  finite  element  method  in  three  dimensions  was  seen  in  section  3.4.2.3  to  have  an 
error  on  the  order  of  2-3  dB  in  the  comparisons  with  the  modal  series  solution  for  a  rigid/fixed 
sphere.  This  suggests  that  the  differences  in  predicted  target  strength  for  some  of  the  side  lobes 
are  not  entirely  an  error  but  possibly  an  effect  of  diffraction.  The  maximum  deviations  in  the 
predictions  generally  increase  with  angle  of  orientation,  suggesting  that  the  effects  of  diffraction 
increase  in  this  region  as  well. 


4.3  Experiment  Background 

There  are  published  data  of  the  forward  scattering  and  backscattering  from  a  sand  dollar 
test  and  a  machined  aluminum  disk  of  similar  size  (Stanton  and  Chu,  2004).  Data  for  scattering 
from  a  bivalve  were  also  collected  in  that  study  but  are  not  described  here.  In  the  experiment, 
the  targets  were  placed  individually  in  a  large  freshwater  tank  and  subjected  to  broadband 
chirps  over  the  frequency  range  40-95  kHz.  The  targets  were  each  rotated  in  the  horizontal 
plane  in  increments  of  one  degree  over  360  degrees  around  a  semi-major  axis  perpendicular  to 
the  direction  of  the  incident  wave.  The  results  were  processed  in  real  time  to  remove  any  multi- 
path  echoes  from  the  walls  of  the  tank.  The  target  information  is  listed  in  table  1. 
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Target 

Horizontal  dimension 
(cm) 

Vertical  dimension 
(cm) 

Upper  camber 
(cm) 

Sand  dollar  ( Dendraster 
excentricus ) 

7.25 

6.7 

.1 

Bivalve  (Dinocar dium 
robustum  vanhyningi ) 

7.0 

6.9 

2.9 

Aluminum  disk 

8.0 

8.0 

0.19 

Table  1  -  Target  information  for  the  acoustic  scattering  experiment  performed  by  Stanton  and  Chu 
(2004). 


The  upper  camber  measurement  in  table  1  refers  to  the  height  of  the  targets  when  placed  on  a 
flat  surface  (sand  dollar  on  its  flat  side  and  bivalve  on  its  concave  side).  The  value  for  the 
aluminum  disk  refers  to  its  thickness,  1.9  mm.  Computed  tomography  (CT)  scans  were  also 
obtained  for  the  sand  dollar  at  the  Marine  Research  Facility  at  the  Woods  Hole  Oceanographic 
Institution,  illustrating  well  the  interior  of  the  test. 


Figure  39  -  Cross  section  of  a  sand  dollar  (Dendraster  excentricus)  obtained  from  CT  scans.  Details  of 
the  inner  structure  are  revealed  in  addition  to  the  high-resolution  measurements  of  the  outer  shape 
(from  Andone  Lavery,  personal  communication). 
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4.4  Comparisons  of  Acoustic  Scattering  Predictions  to 
Experimental  Data 

The  predicted  scattering  based  on  the  Kirchhoff  and  finite  element  methods  can  be 
compared  to  the  experimental  data  obtained  by  Stanton  and  Chu  (2004).  Each  comparison  is 
performed  over  a  range  of  angles  of  orientation  at  a  frequency  of  70  kHz.  First  the  experimental 
results  for  the  aluminum  disk  are  compared  to  the  analytic  solution  of  the  Kirchhoff  integral  for 
a  rigid/ fixed  disk  and  the  finite  element  solution  for  a  rigid/fixed  disk.  Next,  the  experimental 
results  for  scattering  from  the  flat  side  of  the  sand  dollar  are  compared  with  the  analytic 
solution  of  the  Kirchhoff  integral  for  a  rigid / fixed  disk,  the  numerical  solution  of  the  Kirchhoff 
integral  for  the  flat  side  rigid/fixed  sand  dollar  model  (which  uses  the  actual  shape  of  the  sand 
dollar  obtained  from  CT  scans),  and  the  finite  element  method  solution  for  a  rigid/fixed  disk. 
Finally,  the  experimental  results  for  scattering  from  the  round  side  of  the  sand  dollar  are 
compared  to  the  analytic  solution  of  the  Kirchhoff  integral  for  a  rigid/fixed  spherical  cap  and 
the  numerical  solution  of  the  Kirchhoff  integral  for  the  round  side  rigid/fixed  sand  dollar 
model.  Time  limitations  prevented  the  comparison  of  the  finite  element  method  solution  for  the 
rigid/ fixed  sand  dollar  model  (either  flat  side  or  round  side)  or  rigid/ fixed  spherical  cap  with 
experimental  data.  Tables  2,  3,  and  4  summarize  the  comparisons  and  there  location  within  this 
work. 
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Aluminum  Disk 


Kirchhoff, 

analytic 

Kirchhoff, 

numerical 

Finite  element 
method 

Experimental 

Data 

Kirchhoff, 

analytic 

NA 

Section  2.7.31 

Section  4.21 

Section  4.4.1* 

Kirchhoff, 

numerical 

Section  2.7.31 

NA 

None2 

None2 

Finite  element 
method 

Section  4.21 

None2 

NA 

Section  4.4.1* 

Experimental 

Data 

Section  4.4.1* 

None2 

Section  4.4.1* 

NA 

Table  2  -  Comparison  summary  for  aluminum  disk.  Comparisons  were  for  rigid/fixed  disks,  but  with 
different  dimensions  than  the  aluminum  disk.  2The  analytic  and  numerical  solutions  of  the  Kirchhoff 
integral  for  a  rigid/fixed  disk  are  identical,  so  to  avoid  redundancy,  these  comparisons  were  not  made. 
'These  are  revisited  in  section  4.5.2  with  an  added  heuristic  correction,  known  as  the  penetrable 
condition. 


Sand  Dollar,  Flat  Side 


Kirchhoff, 

analytic1 

Kirchhoff, 

numerical 

Finite  element 
method1 

Experimental 

Data 

Kirchhoff, 

analytic1 

NA 

Section  4.4.2.2 

Section  4.2 

Section  4.4.2.2 

Kirchhoff, 

numerical 

Section  4.4.2.2 

NA 

None2 

Section  4.4.2.2* 

Finite  element 
method1 

Section  4.2 

None2 

NA 

Section  4.4.2.2* 

Experimental 

Data 

Section  4.4.2.2 

Section  4.4.2.2* 

Section  4.4.2.2‘ 

NA 

Table  3  -  Comparison  summary  for  sand  dollar,  flat  side.  Models  used  a  rigid/fixed  disk.  2This 
comparison  is  not  shown  due  to  the  similarity  of  the  analytic  and  numerical  solutions  of  the  Kirchhoff 
integral.  'These  are  revisited  in  section  4.5.2  with  an  added  heuristic  correction,  known  as  the 
penetrable  condition. 


103 


Sand  Dollar,  Round  Side 


Kirchhoff, 

analytic1 

Kirchhoff, 

numerical 

Finite  element 
method2 

Experimental 

Data 

Kirchhoff, 

analytic1 

NA 

Section  4.4.23 

None 

Section  4.4.23 

Kirchhoff, 

numerical 

Section  4.4.2.3 

NA 

None 

Section  4.4.23* 

Finite  element 
method2 

None 

None 

NA 

None 

Experimental 

Data 

Section  4.4.23 

Section  4.4.23* 

None 

NA 

Table  4  -  Comparison  summary  for  sand  dollar,  round  side,  ^odel  used  a  rigid/fixed  spherical  cap, 
and  analytic  integral  is  solved  numerically.  2Finite  element  method  simulations  for  a  spherical  cap 
were  not  run  due  to  insufficient  resources.  This  is  revisited  in  section  4.5.2  with  an  added  heuristic 
correction,  known  as  the  penetrable  condition. 


4.4.1  Aluminum  Disk 

The  aluminum  disk  experimental  results  are  first  compared  with  the  analytic  solution  of 
the  Kirchhoff  integral  for  a  rigid/fixed  disk  (equation  (2.48)).  Recall  that  the  numerical  and 
analytic  solutions  of  the  Kirchhoff  integral  for  a  disk  are  identical,  and  so  the  results  for  the 
numerical  Kirchhoff  solution  are  not  shown  here.  The  model  disk  was  given  the  same 
dimensions  as  the  aluminum  disk  from  the  experiment,  namely  a  radius  of  4  cm  and  a  thickness 
of  1.9  mm.  Figure  40  shows  the  Kirchhoff  solution  and  the  experimental  results  for  the 
aluminum  disk  plotted  over  a  range  of  angles  of  orientation. 
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Figure  40  -  Comparison  of  measured  target  strength  at  70  kHz  as  a  function  of  angle  of  orientation  for 
an  aluminum  disk  (solid  starred  line),  radius  4  cm  and  thickness  1.9  mm,  to  predicted  target  strength 
of  a  rigid/fixed  disk  of  similar  dimensions  found  by  solving  the  Kirchhoff  integral  analytically 
(dashed  line). 

Near  broadside  (approximately  +/-  20  degrees),  there  is  good  structural  agreement 
between  the  predictions  and  the  experimental  results.  The  main  lobe,  first  side  lobes,  and  most 
of  the  second  side  lobes  are  visible  in  both  the  analytic  Kirchhoff  solution  and  the  experimental 
results.  However,  the  Kirchhoff  integral  solution  predicted  higher  values  for  the  peaks  of  these 
lobes,  with  the  main  peak  of  the  Kirchhoff  prediction  almost  3.7  dB  above  the  experimental 
results.  The  source  of  this  error  stems  from  either  the  Kirchhoff  method  itself  or  the  experiment. 
Neglecting  experimental  error,  the  two  possible  sources  of  error  in  the  Kirchhoff  method  are  the 
neglect  of  diffraction  and  the  use  of  a  rigid/ fixed  boundary  condition.  Because  of  the  high 
aspect  ratio  shape  near  broadside,  the  effects  of  elastic  waves  are  lessened  at  this  orientation. 


Disk  (Kirchhoff  method  and  experiment) 
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However,  the  rigid/fixed  boundary  condition  still  causes  errors  because  it  neglects  transmission 
of  sound  into  the  object.  In  reality,  not  all  of  the  incident  sound  wave  should  be  reflected  from 
the  front  surface  of  the  disk  as  the  rigid/ fixed  boundary  condition  assumes;  some  of  it  is 
transmitted  through  the  front  surface  and  is  either  reflected  from  the  back  of  the  disk  or 
continues  all  the  way  through  it.  The  discrepancy  seems  to  be  responsible  for  the  Kirchhoff 
method  predicting  the  higher  return  at  broadside.  At  around  25  degrees  away  from  broadside, 
the  peak  null  patterns  lose  their  similarity.  The  effects  of  elastic  waves  and/  or  diffraction  are 
likely  very  important  in  this  region,  as  the  Kirchhoff  method  solution  fails  to  match  the  null 
pattern  of  the  experimental  results. 

Next,  the  experimental  results  for  the  aluminum  disk  are  compared  with  the  predictions 
from  the  finite  element  method  for  a  rigid/ fixed  disk.  Scattering  simulations  for  a  disk  of  radius 
4  cm  and  thickness  1.9  mm  were  run  in  COMSOL,  with  SMD  =  20  and  EPW  =  6  .  The 
frequency  was  set  at  70  kHz  to  match  the  experimental  results  available,  and  the  disk  was 
rotated  in  increments  of  two  degrees  over  a  range  of  angles.  The  far-field  target  strength 
predictions  were  made  with  the  spherical  integral  method  from  Chapter  3.  Figure  41  shows  the 
finite  element  method  prediction  and  the  experimental  scattering  results  for  the  aluminum  disk 
plotted  over  a  range  of  angles  of  orientation. 
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Disk  (finite  element  method  and  experiment) 


Figure  41  -  Comparison  of  measured  target  strength  at  70  kHz  as  a  function  of  angle  of  orientation  for 
an  aluminum  disk  (solid  starred  line),  radius  4  cm  and  thickness  1.9  mm,  to  predicted  target  strength 
of  a  rigid/fixed  disk  of  similar  dimensions  using  the  finite  element  method  in  COMSOL  (thin  solid 
line  with  x's).  EPW  =  6  and  SMD  =  20. 


Like  the  Kirchhoff  method,  the  finite  element  method  predictions  capture  the  main 
features  of  the  experimental  curve  near  broadside.  Also  similar  to  the  Kirchhoff  predictions,  the 
finite  element  method  predicts  larger  target  strength  than  the  experimental  results  in  this  region 
and  is  5.1  dB  greater  at  broadside.  Because  the  finite  element  method  incorporates  the  effects  of 
diffraction,  this  strongly  suggests  that  the  larger  returns  near  broadside  in  its  prediction  are  due 
to  the  use  of  rigid/ fixed  boundary  conditions.  The  finite  element  prediction  is  able  to  capture 
some  of  the  experimental  features  at  higher  angles,  particularly  the  peak  near  sixty  degrees,  and 
to  a  lesser  extent,  the  peak  near  thirty  degrees,  although  they  are  slightly  shifted  in  angle.  The 
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Kirchhoff  method  did  not  accurately  reproduce  the  experimental  data  in  this  region,  suggesting 
that  there  are  strong  diffraction  effects  from  the  disk  at  higher  angles  of  incidence. 


4.4.2  Sand  Dollar 

Unlike  the  aluminum  disk,  the  sand  dollar  does  not  have  an  analytic  solution  to  the 
Kirchhoff  integral.  However,  its  complex  geometry  may  be  approximated  by  simpler  shapes, 
such  as  a  disk  for  the  flat  side  and  a  spherical  cap  for  the  round  side,  which  do  have  analytic 
solutions  to  the  Kirchhoff  integral  (Chapter  2).  Another  approach,  described  in  the  next  section, 
is  to  use  the  numerical  method  of  solving  the  Kirchhoff  integral  by  creating  an  accurate  surface 
mesh  of  the  sand  dollar,  obtained  from  CT  scans.  In  the  two  sections  following,  the  analytic 
solutions  of  the  Kirchhoff  integral  for  the  simple  shapes  are  compared  with  the  numerical 
Kirchhoff  solutions  this  new  sand  dollar  model.  This  was  not  necessary  for  the  disk,  whose 
analytic  and  numerical  solutions  of  the  Kirchhoff  integral  were  shown  to  match  in  section  2.7.3. 
Next  the  experimental  scattering  results  for  the  sand  dollar  are  compared  with  predictions  from 
the  appropriate  models,  described  in  each  section. 
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4.4.2.1  Surface  Mesh  Geometry  from  CT  Scans 

CT  scans  were  performed  on  the  same  sand  dollar  that  was  used  in  the  scattering 
experiment  (Stanton  and  Chu,  2004)  and  were  used  to  create  the  surface  mesh  necessary  for 
solving  the  Kirchhoff  integral  numerically.  Using  the  model  reconstruction  software  AMIRA®, 
the  CT  scans  were  recreated  into  a  three  dimensional  model  of  the  sand  dollar.  This  was 
imported  into  COMSOL  and  the  surface  mesh  information  was  gathered  for  the  numerical 
solution  of  the  Kirchhoff  integral. 


Figure  42  -  Three  dimensional  sand  dollar  model  in  COMSOL,  with  the  round  side  visible. 

In  the  following  two  sections,  the  analytic  solutions  of  the  Kirchhoff  integral  for  the 
simple  shapes  are  first  compared  with  the  numerical  solutions  for  the  sand  dollar.  Then  each  is 
compared  to  the  experimental  data.  This  was  not  necessary  for  the  disk,  whose  numerical  and 
analytic  solutions  of  the  Kirchhoff  integral  were  shown  to  match  in  section  2.7.3. 
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4A.2.2  Sand  Dollar  Flat  Side 


The  experimental  scattering  results  from  the  sand  dollar  are  divided  into  two  categories, 
defined  by  the  sand  dollar  face  oriented  toward  the  incident  wave:  flat  side  and  round  side.  A 
disk  was  used  to  model  the  flat  side  for  an  analytic  solution  to  the  Kirchhoff  integral,  presented 
as  equation  (2.48).  The  radius  of  the  disk  was  3.625  cm,  half  of  the  sand  dollar's  horizontal 
dimension,  and  its  thickness  was  5.5  mm,  half  of  the  maximum  of  the  sand  dollar's  variable 
thickness.  The  choice  of  thickness  was  done  to  create  good  agreement  between  the  analytic  and 
numerical  solutions  of  the  Kirchhoff  integral.  The  Kirchhoff  integral  was  solved  numerically  for 
the  surface  of  the  sand  dollar  itself,  using  the  results  from  the  CT  scans.  Figure  43  shows  a 
coarse  version  of  the  surface  mesh  in  MATLAB,  after  implementing  Newell's  algorithm, 
oriented  with  the  flat  side  up. 


x  10-3 


Figure  43  -  Sand  dollar  surface  mesh  in  MATLAB  after  implementing  Newell's  algorithm,  with  the 
orientation  corresponding  to  an  incident  wave  normal  to  the  flat  side. 
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By  comparing  the  two  methods  of  solving  the  Kirchhoff  integral,  the  effects  of  simplifying 
the  geometry  to  a  disk  in  the  analytic  expression  can  be  seen.  Target  strength  for  these  two 
models  is  shown  in  figure  44  versus  angle  of  orientation,  where  zero  is  normal  to  the  flat  side. 
The  axis  of  rotation  is  the  y  axis  as  seen  in  figure  42. 


Disk  (analytic  Kirchhoff  method) 
and  sand  dollar  flat  side  (numerical  Kirchhoff  method) 


angle  of  orientation 


Figure  44  -  Comparison  of  predicted  target  strength  at  70  kHz  as  a  function  of  angle  of  orientation  for 
a  rigid/fixed  disk,  radius  3.625  cm  and  thickness  5.5  mm,  based  on  the  analytic  solution  of  the 
Kirchhoff  integral  (dashed  line)  to  predicted  target  strength  of  the  flat  side  of  a  rigid/fixed  sand  dollar 
of  similar  dimensions  based  on  a  numerical  solution  to  the  Kirchhoff  integral  (solid  line  with  circles). 


There  is  very  good  agreement  over  a  wide  range  of  angles  for  the  two  methods,  with  deviations 
starting  to  appear  near  forty  degrees  from  broadside.  These  discrepancies  are  due  simply  to 
geometric  differences  between  the  disk  and  the  flat  side  of  the  sand  dollar.  Thus,  a  rigid/  fixed 
disk  proves  to  be  a  good  simplification  of  the  flat  side  of  the  sand  dollar  when  solving  the 
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Kirchhoff  integral.  The  experimental  data  from  the  flat  side  of  the  sand  dollar  is  compared  with 
the  analytic  solution  of  the  Kirchhoff  integral  for  a  rigid/fixed  disk  in  figure  45  and  the 
numerical  solution  of  the  Kirchhoff  integral  for  the  flat  side  of  the  sand  dollar  in  figure  46  over  a 
range  of  angles  of  orientation. 


Disk  (Kirchhoff  method)  and  sand  dollar  flat  side  (experiment) 


angle  of  orientation 


Figure  45  -  Comparison  of  measured  target  strength  at  70  kHz  as  a  function  of  angle  of  orientation  for 
the  flat  side  of  a  sand  dollar  (solid  starred  line),  to  predicted  target  strength  of  a  rigid/fixed  disk, 
radius  3.625  cm  and  thickness  5.5  mm,  found  by  solving  the  Kirchhoff  integral  analytically  (dashed 
line). 
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Figure  46  -  Comparison  of  measured  target  strength  at  70  kHz  as  a  function  of  angle  of  orientation  for 
the  flat  side  of  a  sand  dollar  (solid  starred  line),  to  predicted  target  strength  of  the  flat  side  of  a 
rigid/fixed  sand  dollar  found  by  solving  the  Kirchhoff  integral  numerically  (solid  line  with  circles). 

Because  the  two  methods  of  solving  the  Kirchhoff  integral  produced  such  similar  results, 
the  focus  will  stay  on  figure  46.  While  the  lobe  structures  of  the  Kirchhoff  and  experimental 
curves  are  similar  with  respect  to  angle  of  orientation,  there  are  significant  differences  between 
their  target  strength  values  at  the  locations  of  the  peaks  and  nulls.  Generally,  the  experimental 
data  is  below  the  predicted  results  from  the  Kirchhoff  method.  The  main  lobe  of  the  experiment 
is  8.1  dB  below  the  main  lobe  of  the  Kirchhoff  curve.  This  is  a  similar  result  to  the  comparison 
with  the  aluminum  disk.  Most  likely,  modeling  the  sand  dollar  as  a  rigid/fixed  scatterer  is  the 
cause  of  the  higher  predicted  target  strength,  for  the  reasons  mentioned  before.  It  is  important 
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to  note,  however,  the  Kirchhoff  solution  maintains  lobe  structure  similarity  with  the  aluminum 
disk  over  the  entire  range  of  angles,  unlike  the  aluminum  disk. 

Next,  the  experimental  results  for  the  flat  side  of  the  sand  dollar  are  compared  with  the 
predictions  from  the  finite  element  method,  in  which  the  flat  side  of  the  sand  dollar  was 
modeled  as  a  rigid/fixed  disk.  Scattering  simulations  from  a  disk  with  radius  3.625  cm  and 
thickness  5.5  mm  were  run  in  COMSOL  with  SMD  =  20  and  EPW  =  6  .  The  simulations  were 
performed  at  a  frequency  of  70  kHz  and  as  a  function  of  orientation,  in  which  the  disk  was 
rotated  in  increments  of  one  degree  over  a  range  of  angles.  The  far-field  calculations  were 
obtained  from  the  spherical  integral  method  from  Chapter  3.  Figure  47  shows  the  finite  element 
solution  and  the  experimental  scattering  results  from  the  sand  dollar's  flat  side  over  a  range  of 
angles. 
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Disk  (finite  element  method)  and  sand  dollar  flat  side  (experiment) 


Figure  47  -  Comparison  of  measured  target  strength  at  70  kHz  as  a  function  of  angle  of  orientation  for 
the  flat  side  of  a  sand  dollar  (solid  starred  line),  to  predicted  target  strength  of  a  rigid/fixed  disk  of 
radius  3.625  cm  and  thickness  5.5  mm  using  the  finite  element  method  in  COMSOL  (thin  solid  line 
with  x's).  EPW  =  6  and  SMD  =  20. 


The  predictions  from  the  finite  element  method  are  again  greater  than  the  experimental 
results  at  low  angles  of  incidence.  At  broadside,  the  finite  element  curve  is  about  9.2  dB  greater 
than  the  experimental  value.  Again,  there  is  not  much  difference  in  the  angular  dependence  of 
the  peak  and  nulls  between  the  two  curves. 


4.4.2.3  Sand  Dollar  Round  Side 

A  spherical  cap  was  used  for  the  solution  of  the  Kirchhoff  integral  for  the  round  side  of 
the  sand  dollar.  To  create  the  cap,  a  sphere  of  radius  6.523  cm  defined  for  6  =  0  to  6  =  33.8 


115 


degrees  was  used.  This  gave  a  cap  with  radius  3.625  cm  and  thickness  1.1  cm,  the  dimensions  of 
the  sand  dollar.  The  expression  for  a  spherical  cap,  equation  (2.50),  was  evaluated  with  the 
midpoint  rule  in  MATLAB.  This  process  must  not  be  confused  with  the  numerical  solution  of 
the  Kirchhoff  integral  using  the  sand  dollar  surface  mesh.  For  this  numerical  solution,  the  round 
side  of  the  sand  dollar  was  modeled  again  using  the  CT  scans,  with  a  coarse  mesh  of  the  round 
side  seen  in  figure  48. 


Figure  48  -  Sand  dollar  surface  mesh  in  MATLAB  after  implementing  Newell's  algorithm,  with  the 
orientation  corresponding  to  an  incident  wave  normal  to  the  round  side. 

Again,  these  two  methods  of  solving  the  Kirchhoff  integral  were  compared  to  see  the 
validity  of  the  geometry  simplification.  The  target  strengths  of  the  two  solutions  are  shown  in 
figure  49  over  a  range  of  angles  of  orientation,  where  zero  corresponds  to  broadside.  The  sand 
dollar  model  was  again  rotated  around  the  y  axis  as  illustrated  in  figure  42. 
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Spherical  cap  (Kirchhoff  method) 
and  sand  dollar  round  side  (numerical  Kirchhoff  method) 
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Figure  49  -  Comparison  of  predicted  target  strength  at  70  kHz  as  a  function  of  angle  of  orientation  for 
an  spherical  cap,  radius  3.625  cm  and  thickness  1.1  cm,  based  on  the  analytic  solution  of  the  Kirchhoff 
integral  (dashed  line)  to  predicted  target  strength  of  the  round  side  of  a  rigid/fixed  sand  dollar  of 
similar  dimensions  based  on  a  numerical  solution  to  the  Kirchhoff  integral  (solid  line  with  circles). 


The  two  predictions  have  some  general  similarities  but  also  some  very  distinct 
differences.  The  agreement  between  the  two  methods  is  worse  than  that  between  the  sand 
dollar's  flat  side  and  the  disk.  The  discrepancies  are  most  noticeable  at  angles  close  to 
broadside,  where  the  spherical  cap  approach  predicted  a  null,  but  the  numerical  approach 
predicted  a  peak.  Interestingly,  there  is  far  less  structure  at  angles  far  from  broadside  for  the 
models  of  the  sand  dollar's  round  side  as  opposed  to  those  for  its  flat  side.  However,  at  angles 
larger  than  approximately  15  degrees,  the  agreement  between  these  two  models  was  relatively 
good,  with  the  target  strength  decreasing  with  angle  and  the  levels  in  relatively  good 
agreement. 
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The  experimental  scattering  data  from  the  sand  dollar's  round  side  are  compared  with  the 
analytic  solution  of  the  Kirchhoff  method  for  the  spherical  cap  in  figure  50  and  the  numerical 
solution  for  the  sand  dollar  model  in  figure  51. 


Spherical  cap  (Kirchhoff  method)  and  sand  dollar  round  side  (experiment) 


Figure  50  -  Comparison  of  measured  target  strength  at  70  kHz  as  a  function  of  angle  of  orientation  for 
the  round  side  of  a  sand  dollar  (solid  starred  line),  to  predicted  target  strength  of  a  rigid/fixed 
spherical  cap,  radius  3.625  cm  and  thickness  1.1  cm,  found  by  solving  the  Kirchhoff  integral 
analytically  (dashed  line). 
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Sand  dollar  round  side  (Kirchhoff  method  and  experiment) 


Figure  51  -  Comparison  of  measured  target  strength  at  70  kHz  as  a  function  of  angle  of  orientation  for 
the  round  side  of  a  sand  dollar  (solid  starred  line),  to  predicted  target  strength  of  the  round  side  of  a 
rigid/fixed  sand  dollar  found  by  solving  the  Kirchhoff  integral  numerically  (solid  line  with  circles). 


Again,  the  predictions  based  on  the  analytic  (figure  50)  and  numerical  (figure  51) 
solutions  of  the  Kirchhoff  integral  mimic  the  experimental  results  in  their  general  structure  but 
have  noticeable  and  important  differences.  The  central  peak  at  broadside  incidence  observed  in 
the  experimental  data  supports  the  hypothesis  that  the  numerical  Kirchhoff  solution  is  a  better 
model  than  the  spherical  cap.  Once  again,  the  Kirchhoff  results  have  higher  values  than  the 
experimental  data,  likely  caused  by  the  choice  of  rigid/ fixed  boundary  conditions.  While  it  is 
probable  that  some  fraction  of  the  incident  wave  is  transmitted  through  the  sand  dollar  in  the 
actual  experiment,  this  is  not  accounted  for  with  the  rigid/ fixed  boundary  condition.  The 
increased  structure  in  the  lobes  of  the  experiment  is  likely  caused  by  diffraction  effects. 
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Unfortunately,  due  to  time  constraints  and  long  numerical  run  times,  the  finite  element  method 
was  not  tested  for  a  spherical  cap,  making  this  hypothesis  difficult  to  verify. 

4.5  Heuristic  Improvement  to  Scattering  Models 
Developed  for  Rigid/ Fixed  Objects 

There  are  errors  inherent  in  modeling  an  elastic  scatterer  as  a  rigid/  fixed  object  because 
the  portion  of  the  wave  that  is  transmitted  through  the  front  surface  of  the  scatterer  is  neglected. 
In  addition  to  neglecting  the  transmitted  portion  of  the  incident  wave,  the  models  presented  in 
this  thesis  for  aluminum  disks  and  sand  dollars  have  neglected  the  effects  of  elastic  shear 
waves.  The  effects  of  elastic  shear  waves  are  extremely  complicated  for  irregular  shapes  such  as 
the  sand  dollar  and  are  beyond  the  scope  of  this  work.  Instead,  the  focus  of  the  next  section  will 
be  on  attempting  to  quantify  how  much  of  the  incident  wave  is  transmitted  through  versus 
reflected  from  the  sand  dollar.  This  proves  easiest  when  the  high  aspect  ratio  shape  is  near 
broadside.  As  will  be  shown,  a  heuristic  correction  to  the  rigid/ fixed  boundary  condition  used 
for  modeling  the  aluminum  disk  and  sand  dollar  can  be  made  that  works  very  well  near 
broadside  incidence. 

4.5.1  Reflection  Coefficients  from  Infinite  Half-Spaces  and  Layers 

When  a  wave  arrives  at  a  boundary,  some  of  it  will  be  reflected  and  some  will  be 
transmitted  into  the  new  medium.  If  the  boundary  is  planar,  the  reflection  coefficient  can 
quantify  how  much  of  the  incident  wave  is  reflected.  If  an  incident  wave  has  amplitude  pQ ,  the 
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reflected  wave  has  amplitude  Rp0 .  In  the  case  of  a  rigid/ fixed  boundary,  the  entire  wave  is 

reflected,  none  is  transmitted,  and  R  =  1 .  A  simple  case  for  determining  R  is  from  the  planar 
boundary  of  two  fluid  infinite  half  spaces.  When  a  plane  wave  traveling  through  medium  0 
(c0 ,  p0 )  arrives  at  angle  0  to  the  normal  of  the  boundary  with  medium  1  (c,  ,px),  the  reflection 
coefficient  is  defined  as 


R 


m  cos 


0-yJn2  -  sin2  0 


(4.1) 


m cos  0  +  V«2  -  sin2  6 


where  n  =  c0  / 1\  and  m  =  pj  p0  (Frisk,  1994).  At  normal  incidence,  this  becomes 


m  -  n  (4.2) 

R  — - . 

m  +  n 

Now  consider  the  reflection  off  of  a  fluid  layer  in  between  two  fluid  infinite  half  spaces  at 
normal  incidence.  A  wave  travels  through  medium  0  (c0,p0)  and  arrives  perpendicular  to  the 
boundary  with  medium  1  (q ,  px ) .  Medium  1  forms  a  layer  of  thickness  h ,  beyond  which  is 
medium  2  (c2,p2) . 
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Figure  52  -  Geometry  of  a  single  layer  in  between  two  infinite  half  spaces. 

There  are  an  infinite  number  of  contributions  to  the  total  reflection  from  this  layer:  the  wave 
that  reflects  from  the  front  of  the  layer  and  waves  that  reflect  within  the  layer  any  number  of 
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times  before  transmitting  back  into  medium  0.  Brekhovskikh  and  Lysanov  (2003)  give  a  closed 
form  expression  for  the  layer's  reflection  coefficient. 


Rm+Rl2e^h  ^-3) 

L  1  +  RmRve2ik'h' 

where  Rab  is  the  reflection  coefficient  of  a  wave  in  medium  a  off  of  medium  b  .  The 
wavenumber  within  the  layer  is  kx .  A  similar  calculation  can  be  done  if  the  fluid  layer  is 

replaced  by  an  infinite  elastic  plate.  If  the  incident  wave  arrives  normal  to  the  plate,  then  there 
is  no  shear  stress  and  only  compressional  waves  arise.  Then  the  elastic  material  may  be  treated 
as  a  fluid. 

Now  imagine  an  incident  plane  wave  arriving  normal  to  n  layers.  The  j  th  layer  has 
thickness  hj  and  material  properties  c  .  and  p  ■ .  The  total  reflection  coefficient  from  the  series 

of  layers  can  be  found  using  a  recursive  process.  The  total  reflection  from  layer  n  is  found  first 
using  a  modified  version  of  the  expression  from  Brekhovskikh  and  Lysanov  (2003), 

R  .  +R  +.e2iKK  (4-4) 

f)  _  n—l,n  n,n+ 1 

L-n  I  i  O  p  „2* A 

1  +  Kn-l,nKn,n+ie 

This  can  then  be  used  to  find  the  total  reflection  from  both  layers  n  and  n—  1 , 

Rn  2(,  ,+RLne2lk"-'h"-'  (4-5) 

r>  _  n— Z,n— 1  L,n 

i’”“1  l  +  R  R  e^K-A-i  ■ 

The  process  can  be  continued  for  each  layer  until  the  total  reflection  is  calculated.  As  was  the 
case  with  one  layer,  at  normal  incidence,  this  process  can  be  used  for  fluid  or  elastic  layers,  or  a 


combination  of  the  two. 
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4.5.2  Determination  of  Reflection  Coefficient  for  Aluminum  Disk  and 
Sand  Dollar 


Now  the  scattering  from  a  finite  bounded  object  is  heuristically  modeled  by  assigning  it 
a  reflection  coefficient  found  in  the  same  manner  as  the  reflection  coefficients  of  fluid  and  elastic 


layers.  The  scattering  from  an  object  is  approximated  as 


ikr 


P scat  RLP0  fimpen^Q')  ’ 


(4.6) 


where  RL  is  its  reflection  coefficient  and  fh  is  the  scattering  amplitude  of  the  object  if  it  had  a 

rigid/ fixed  boundary.  This  approximation  assumes  the  reflection  coefficient  is  constant 
everywhere  on  the  object,  which  is  inaccurate  for  many  shapes  due  to  different  angles  of 
incidence.  However,  this  assumption  is  reasonable  at  broadside  for  a  high  aspect  ratio  shape 
such  as  a  sand  dollar  or  disk.  In  the  backscattering  direction  with  this  formulation, 

rs  =  ioiog|/;j;J  (47) 

TS  =  101og|/.ffl(leJ(|2  +  101og|/?J' . 


The  approximation  gives  that  the  target  strength  of  an  object  is  equal  to  its  rigid/ fixed  target 
strength  plus  a  term  that  is  a  function  of  its  reflection  coefficient.  The  rigid/ fixed  target  strength 
is  what  was  calculated  using  both  the  Kirchhoff  and  finite  element  methods.  These  results  can 
easily  be  manipulated  to  include  the  effects  of  the  reflection  coefficient. 

Approximations  for  reflection  coefficients  of  the  sand  dollar  and  aluminum  disk  can  be 
made  by  finding  them  for  infinite  plates  with  similar  thicknesses  and  material  properties.  This 
technique  should  work  well  near  broadside,  where  the  dimensions  of  the  objects  allow  such  an 
approximation.  Away  from  broadside,  this  method  is  not  expected  to  perform  well.  The 
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aluminum  disk  can  be  modeled  as  an  infinite  plate  of  thickness  1.9  mm  with  water  ( p  =  1000 


kg/ m3,  c  =  1500  m/s)  on  either  side.  The  values  for  the  compressional  sound  speed  and  density 
for  aluminum  were  taken  from  Molz  and  Beamish  (1996). 


Material 

Thickness 

(mm) 

Sound  speed 
(compress.) 

Density 

Rl 

|Rl| 

lOlog  |  Rl  |  2 

Aluminum 

1.9  mm 

5360  m/  s 

2700  kg/  m3 

0.362  - 
0.472i 

0.5954 

-4.50  dB 

Table  5  -  Reflection  coefficient  calculation  for  aluminum  disk. 


This  suggests  that  modeling  the  aluminum  disk  as  a  rigid/fixed  scatterer  added  more  than  4  dB 
to  its  true  target  strength  at  broadside. 

The  sand  dollar  is  more  complicated  to  model  because  of  its  irregular  shape.  The  CT 
scans  of  the  sand  dollar  (e.g.  figure  39)  provide  clues  on  how  to  simplify  it.  The  top  and  bottom 
interfaces  are  very  distinct,  separately  by  a  fluid  interior.  Ignoring  both  the  curvature  of  the  top 
and  the  outer  edge  regions  allows  the  sand  dollar  to  be  modeled  as  two  infinite  plates,  with 
fluid  layer  separating  them.  The  thicknesses  of  these  three  layers  must  be  estimated  from  the  CT 
scan  itself.  The  maximum  camber  of  the  sand  dollar  is  1.1  cm,  so  this  places  an  upper  bound  on 
the  width  of  the  three  layers.  The  sand  dollar's  test  wall  thickness  is  approximately  1  mm  on 
both  the  top  and  bottom  interfaces.  The  sand  dollar  is  composed  of  calcite,  a  polymorph  of 
calcium  carbonate  (CaCOa).  Carmichael  (1982)  gives  the  density  of  calcite  as  2710  kg/  m3  and  the 
compressional  sound  speed  as  6530  m/ s.  The  reflection  coefficient  can  be  calculated  over  a 
range  of  both  fluid  and  test  layer  thicknesses.  Water  ( p  -  1000  kg/ m3,  c  =  1500  m/ s)  was  used 
for  both  the  fluid  layer  and  the  fluid  half  spaces  on  either  side  of  the  three  layers. 
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Reflection  coefficient  for  three  layer  model  (calcite-water-calcite) 
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Figure  53  -  Reflection  coefficients  from  three  layer  model  (calcite-water-calcite).  Calcite  thickness 
refers  to  both  top  and  bottom  layer  thickness. 

Because  the  fluid  layer  thickness  in  the  actual  sand  dollar  varies  over  the  range  plotted 
in  figure  53,  there  is  a  large  range  of  reasonable  reflection  coefficients,  approximately  0.2  to  0.7. 
In  order  to  choose  the  proper  value  for  \RL\,  the  rigid/fixed  predicted  target  strengths  from  the 
models  are  compared  to  the  experimental  data  near  broadside.  For  the  sand  dollar,  a  reflection 
coefficient  I  RL  1=  0.4  was  chosen  to  match  the  Kirchhoff  prediction  with  the  experimental  data 
for  the  flat  side  at  broadside.  Based  on  figure  53,  this  value  is  reasonable. 
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4.5.3  Heuristic  Correction  to  the  Rigid/Fixed  Scattering  Models 


By  using  the  calculated  reflection  coefficient  for  an  object,  the  heuristically  corrected 
model  values  can  be  calculated  with  equation  (4.7).  The  inclusion  of  the  reflection  coefficient 
into  the  calculations  is  referred  to  here  as  the  penetrable  condition. 

4.5.3.1  Aluminum  Disk 

Using  the  value  of  I RL  1=  0.5954  derived  in  the  previous  section  for  the  reflection 

coefficient  for  the  aluminum  disk,  the  target  strength  based  on  both  the  analytic  solution  to  the 
Kirchhoff  integral  (figure  54)  and  the  finite  element  method  (figure  55)  are  compared  with  the 
experimental  data. 


Disk  (analytic  Kirchhoff  method  with  penetrable  condition  and  experiment) 


angle  of  orientation 


Figure  54  -  Comparison  of  measured  target  strength  at  70  kHz  as  a  function  of  angle  of  orientation  for 
an  aluminum  disk  (solid  starred  line),  radius  4  cm  and  thickness  1.9  mm,  to  predicted  target  strength 
of  a  penetrable  disk  ( |  Rl  |  =  0.5954)  of  similar  dimensions  found  by  solving  the  Kirchhoff  integral 
analytically  (dashed  line). 
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Disk  (finite  element  method  with  penetrable  condition  and  experiment) 


Figure  55  -  Comparison  of  measured  target  strength  at  70  kHz  as  a  function  of  angle  of  orientation  for 
an  aluminum  disk  (solid  starred  line),  radius  4  cm  and  thickness  1.9  mm,  to  predicted  target  strength 
of  an  penetrable  disk  ( |  Rl  |  =  0.5954)  of  similar  dimensions  using  the  finite  element  method  in 
COMSOL  (thin  solid  line  with  x's).  EPW  =  6  and  SMD  =  20. 


With  the  heuristic  correction  that  accounts,  in  part,  for  the  assumption  of  a  rigid/ fixed 
boundary  condition,  the  Kirchhoff  and  finite  element  methods  now  predict  the  aluminum  disk's 
scattering  relatively  well  near  broadside,  especially  along  the  main  lobe  and  first  pair  of  side 
lobes.  The  Kirchhoff  prediction  is  within  1  dB  of  the  experimental  results  at  broadside.  The 
finite  element  method  is  seen  to  have  done  very  well  away  from  broadside  as  well,  capturing 
the  large  side  lobe  at  sixty  degrees. 
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4.5.3.2  Sand  Dollar  Flat  Side 


The  sand  dollar  was  a  more  difficult  shape  for  which  to  determine  reflection  coefficient 
because  of  its  multiple  interfaces  and  differing  thickness.  For  the  flat  side,  a  value  of  I RL  1=  0.4 
was  chosen  and  was  shown  to  be  reasonable  using  the  three  layer  approach.  This  value  was 
chosen  so  the  value  of  the  numerical  solution  to  the  Kirchhoff  integral  at  broadside  would 
closely  match  that  of  the  experimental  data.  Figure  56  shows  the  Kirchhoff  method  prediction 
for  the  sand  dollar's  flat  side  with  the  penetrable  condition,  compared  with  the  experimental 
data  for  the  sand  dollar's  flat  side. 


Sand  dollar  flat  side  (Kirchhoff  method  with  penetrable  condition  and  experiment) 
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Figure  56  -  Comparison  of  measured  target  strength  at  70  kHz  as  a  function  of  angle  of  orientation  for 
the  flat  side  of  a  sand  dollar  (solid  starred  line),  to  predicted  target  strength  of  the  flat  side  of 
penetrable  sand  dollar  ( |  Rl  |  =  0.4)  found  by  solving  the  Kirchhoff  integral  numerically  (solid  line 
with  circles). 
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Figure  57  shows  the  finite  element  method  solution  for  a  disk  with  the  penetrable  condition, 
compared  with  the  experimental  data  for  the  sand  dollar's  flat  side. 


Disk  (finite  element  method  with  penetrable  condition) 


and  sand  dollar  flat  side  (experiment) 


Figure  57  -  Comparison  of  measured  target  strength  at  70  kHz  as  a  function  of  angle  of  orientation  for 
the  flat  side  of  a  sand  dollar  (solid  starred  line),  to  predicted  target  strength  of  a  penetrable  disk  ( |  Rl  | 
=  0.4)  of  radius  3.625  cm  and  thickness  5.5  mm  using  the  finite  element  method  in  COMSOL  (thin 
solid  line  with  x's).  EPW  =  6  and  SMD  =  20. 


Both  model  predictions  are  greatly  improved  with  the  heuristic  correction  that  accounts, 
in  part,  for  the  rigid/fixed  boundary  conditions.  The  Kirch hoff  prediction  is  very  good  near 
broadside  now  with  the  choice  of  reflection  coefficient.  The  finite  element  method  manages  not 
only  to  capture  the  main  lobe  but  also  the  first  pair  of  side  lobes.  The  finite  element  method  is 
higher  than  the  corrected  experimental  curve  because  of  the  reflection  coefficient  chosen.  The 
slight  difference  in  width  of  the  main  lobes  in  figure  57  is  likely  a  result  of  the  choice  of  radius 
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for  the  disk  in  the  finite  element  method  solution.  Both  predictions,  especially  the  finite  element 
solution,  agree  with  the  data  at  higher  angles  of  orientation  also. 


4.5.3.3  Sand  Dollar  Round  Side 

For  the  round  side  of  the  sand  dollar,  the  value  of  the  reflection  coefficient  I RL  1=  0.4 
was  chosen  in  order  to  match  that  of  the  flat  side.  Figure  58  shows  the  numerical  solution  of  the 
Kirchhoff  integral  for  the  sand  dollar  round  side  with  the  penetrable  condition  and  the 
experimental  scattering  results  from  the  sand  dollar  round  side. 


Sand  dollar  round  side 


(numerical  Kirchhoff  method  with  penetrable  condtion  and  experiment) 


angle  of  orientation 


Figure  58  -  Comparison  of  measured  target  strength  at  70  kHz  as  a  function  of  angle  of  orientation  for 
the  round  side  of  a  sand  dollar  (solid  starred  line),  to  predicted  target  strength  of  the  round  side  of  a 
penetrable  sand  dollar  ( |  Rl  |  =  0.4)  found  by  solving  the  Kirchhoff  integral  numerically  (solid  line 
with  circles). 
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There  is  obvious  structural  agreement  between  the  predicted  results  and  the  experimental  data. 
The  very  slight  traces  of  the  nulls  in  the  Kirchhoff  prediction  line  up  with  those  in  the 
experimental  data.  There  is  still  much  less  structure  in  the  model  prediction. 


4.6  Summary  of  Results 

The  scattering  models  developed  in  this  thesis  had  varying  amounts  of  success  in 
predicting  the  experimental  data  for  scattering  from  aluminum  disks  and  sand  dollars.  The 
heuristic  correction,  which  accounted  for  transmission  of  the  incident  wave,  was  developed  and 
was  found  to  significantly  improve  these  predictions.  In  the  case  of  the  aluminum  disk,  both  the 
Kirchhoff  and  finite  element  methods  succeeded  in  accurately  recreating  the  scattering  curves 
near  broadside.  The  finite  element  solution  was  better  away  from  broadside,  suggesting  that 
there  was  strong  diffraction  around  the  disk. 

In  the  case  of  the  flat  side  of  the  sand  dollar,  the  numerical  solution  of  the  Kirchhoff 
integral  proved  to  be  very  similar  to  the  analytic  solution  for  a  disk.  Both  the  Kirchhoff  and 
finite  element  methods  did  a  reasonably  good  job  in  predicting  the  peak-null  structure  of  the 
experimental  data  over  the  entire  range  of  angles.  The  success  of  the  heuristic  correction  over 
the  entire  range  of  orientations  suggests  that  the  sand  dollar  may  potentially  be  modeled  as  a 
fluid  scatterer  in  the  future,  without  the  need  for  elastic  considerations.  Testing  this  is  an 
important  next  step  toward  understanding  the  sand  dollar's  scattering  mechanics.  The  fact  that 
the  finite  element  method  did  not  drastically  outperform  the  Kirchhoff  method  away  from 
broadside,  as  it  did  for  the  aluminum  disk,  suggests  that  the  effects  of  diffraction  were  less  for 
the  sand  dollar;  however,  this  is  difficult  to  determine.  If  it  is  true  that  diffraction  is  not  as 
strong  for  the  sand  dollar  as  the  aluminum  disks,  it  would  be  a  result  of  the  geometries  and/  or 
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the  material  properties  of  the  two  objects.  Understanding  this  has  implications  for  the 
application  of  using  diffraction  of  benthic  shells  to  detect  buried  objects  at  subcritical  angles. 

The  round  side  of  the  sand  dollar  was  not  modeled  as  well  as  the  flat  side  with  the 
Kirchhoff  method.  Because  the  finite  element  method  was  not  tested  for  this  shape,  it  is  difficult 
to  conclusively  determine  the  reasons  for  this.  However,  it  was  seen  that  the  numerical  solution 
to  the  Kirchhoff  integral  is  a  better  model  than  the  analytic  solution  for  a  spherical  cap.  This 
proves  the  utility  of  the  solving  the  Kirchhoff  integral  numerically.  Although  the  results  were 
not  as  good  as  for  the  flat  side,  the  general  structures  of  the  prediction  and  the  experimental 
data  are  similar. 
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Chapter  5  Summary  and  Conclusion 


This  focus  of  this  thesis  is  to  understand  the  acoustic  backscattering  from  individual  sand 
dollars  ( Dendraster  excentricus)  in  the  free  field.  The  following  paragraphs  summarize  the 
motivations  for  this  work,  the  research  performed,  recommendations  for  future  work,  and  the 
contributions  of  this  thesis. 


5.1  Motivation 

Knowledge  of  seafloor  scattering  has  many  important  applications.  Benthic  shelled 
organisms  have  been  shown  to  be  important  contributors  to  scattering  from  the  seafloor  in 
certain  cases,  depending  on  the  acoustic  frequency,  angle  of  incidence,  and  density  of 
organisms.  Understanding  the  scattering  physics  of  individual  benthic  shells  is  an  important 
step  toward  understanding  the  effects  of  collections  in  the  ocean.  Sand  dollars  were  chosen  for 
this  thesis  because  of  existing  laboratory  data  to  which  the  models  could  be  compared. 


5.2  Modeling 

Acoustic  scattering  models  can  provide  insight  to  the  scattering  mechanisms  from 
complex  objects,  such  as  sand  dollars.  The  Kirchhoff  method  is  a  very  simple  model  that  proves 
to  be  very  useful  at  high  frequencies  (short  wavelengths  relative  to  the  object  size).  The 
numerical  method  of  solving  the  Kirchhoff  integral  developed  here  greatly  broadens  this 
approach  by  allowing  the  scattering  from  irregularly  shaped  objects  to  be  predicted.  It  is 
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believed  that  the  accuracy  of  predictions  for  an  object's  scattering  can  be  greatly  increased  by 
solving  the  integral  numerically,  as  opposed  to  solving  it  analytically  for  a  simplified  object,  as 
was  the  sand  dollar's  round  side.  Applied  to  acoustic  scattering  problems,  the  finite  element 
method  is  a  technique  far  more  complicated  than  the  Kirchhoff  method.  It  requires  powerful 
computer  resources  for  any  meaningful  three  dimensional  calculations.  Its  advantage  is  that  it 
does  not  make  any  of  the  approximations  implicit  in  the  Kirchhoff  method,  such  as  neglecting 
the  effects  of  diffraction.  However,  it  has  been  seen  that  at  high  enough  frequencies  near 
broadside,  the  finite  element  method's  predictions  are  of  about  the  same  accuracy  as  those  from 
the  Kirchhoff  method.  It  is  believed  that  the  finite  element  method's  value  in  acoustic  scattering 
problems  lies  at  lower  frequencies.  The  heuristic  correction  for  the  rigid/fixed  boundary 
condition  that  was  developed  here  was  shown  to  improve  the  ability  of  the  model  to  predict  the 
experimental  scattering  data. 


5.3  Recommendations  for  Future  Work 

One  main  limitation  of  this  work  was  the  use  of  a  rigid/ fixed  boundary  condition  in  the 
scattering  models.  This  allowed  only  "first  order"  approximations  to  the  scattering  from  the 
elastic  objects  to  be  studied.  Without  incorporating  the  elasticity  of  a  target,  it  is  difficult  to 
distinguish  the  importance  of  certain  scattering  phenomena  such  as  diffraction,  surface  elastic 
waves,  and  shear  wave  effects.  Although  the  Kirchhoff  method  is  unable  to  incorporate  the 
effects  of  elasticity,  the  same  is  not  true  for  the  finite  element  method.  Currently,  a  finite 
element  code  is  being  developed  at  the  Woods  Hole  Oceanographic  Institution  which  would 
allow  for  the  modeling  of  elastic  domains  (Gonzalo  Feijoo,  personal  communication).  This  code 
will  also  incorporate  a  perfectly  matched  layer,  a  far  better  radiation  condition  than  the 
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simplified  Sommerfeld  condition  used  in  the  COMSOL  simulations  (Berenger  1994, 1996;  Harari 
et  al.,  2000).  It  would  be  valuable  to  test  scattering  predictions  for  sand  dollars  or  other  benthic 
shells  with  the  ability  to  model  their  elasticity.  In  addition,  this  study  concerned  sand  dollars  in 
the  free  field.  Scattering  involving  the  geometry  where  sand  dollars  are  on  or  near  a  boundary 
should  be  examined.  As  sand  dollars  will  generally  be  encountered  in  the  ocean  lying  on  the 
seafloor,  modeling  their  scattering  while  near  a  boundary  could  help  gain  additional  insight. 
This  could  be  performed  in  the  finite  element  method  by  dividing  the  surrounding  medium 
domain  into  two  fluid  regions,  one  representing  the  seawater  and  one  representing  the  seafloor. 

The  comparisons  between  model  predictions  and  experimental  data  in  this  thesis  focused 
on  a  single  frequency  over  a  range  of  angles.  The  insight  that  these  comparisons  offer  could  be 
bolstered  by  further  comparisons  with  more  experimental  data.  Experimental  data  for 
additional  sand  dollars  would  reduce  the  effects  of  experimental  error  as  well  as  provide  insight 
to  their  scattering  physics  by  examining  differences  between  samples.  Comparisons  to  model 
predictions  at  frequencies  other  than  70  kHz  would  aid  this  understanding  as  well.  As  a  final 
experimental  recommendation,  the  scattering  from  a  live  sand  dollar  should  be  studied.  It  is 
unknown  how  the  scattering  physics  of  such  an  organism,  surrounded  by  a  thin  layer  of  tissue, 
would  be  different  from  only  the  test. 

Sand  dollars  are  often  found  in  dense  aggregates  in  situ,  creating  a  very  complicated 
scattering  problem.  The  hope  is  that  the  study  of  scattering  from  individual  sand  dollars  in  the 
free  field  leads  to  a  greater  understanding  of  how  these  aggregates  of  organisms  scatter  sound 
in  the  ocean. 
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5.4  Contributions  of  Thesis 

•  Development  of  analytic  Kirchhoff  method  scattering  models  for  a  sand  dollar,  namely  a 
disk  and  spherical  cap. 

•  Development  of  a  numerical  technique  for  solving  the  Kirchhoff  integral  for  complex 
three  dimensional  objects. 

•  Comparison  of  the  analytic  and  numerical  solutions  of  the  Kirchhoff  integral  with 
experimental  scattering  data  from  a  sand  dollar  and  aluminum  disk. 

•  Application  of  the  finite  element  method  in  the  program  COMSOL  for  predicting  the 
scattering  from  a  rigid / fixed  disk.  Comparison  of  finite  element  method  model 
predictions  with  experimental  scattering  data  from  a  sand  dollar  and  aluminum  disk. 

•  Development  of  a  heuristic  approach  for  accounting  for  a  penetrable  boundary 
condition  for  high  aspect  ratio  objects  at  broadside.  Application  of  this  approach  for  the 
Kirchhoff  and  finite  element  method  predictions  for  scattering  from  a  sand  dollar  and 
aluminum  disk. 
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Appendix  A:  MATLAB  Code 


One  way  to  solve  the  Kirchhoff  integral  numerically  is  by  generating  a  mesh  in 
COMSOL  Multiphysics.  The  object  must  be  oriented  so  that  the  incident  plane  wave  is  traveling 
in  the  negative  z  direction.  Once  the  geometry  and  the  mesh  have  been  created,  they  can  be 
exported  to  COMSOL  script  as  a  fern  structure.  The  function  kirchhoff_mesh.m  was  written  for 
COMSOL  Multiphysics  Script  (ver.  1.0a)  to  extract  the  surface  mesh  information  and  save  to  file 
its  connectivity  and  coordinate  information. 


oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 


;%%  kirchhof f_mesh . m 
;%%  Greg  Dietzen  5/27/08 


oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 


%%%  Creates  a  surface  mesh  of  an  object  from  COMSOL  fern  structure.  Object 
%%%  geometry  and  mesh  must  be  created  in  COMSOL  and  exported  to  COMSOL 
%%%  script  as  an  fern  structure  (Ctrl  f  in  GUI) . 

function  kirchhof f_mesh (fern) ; 

fem.xmesh  =  meshextend ( f em) ; 

%%%  Find  coordinates  of  all  nodes 
nodeinfo  =  xmeshinf o ( fem, ' out nodes ') ; 
nodes  =  [ 0 : length (nodeinfo . dofs ) -1 ] 
coords  =  nodeinfo . coords ' ; 
allnodecoords  =  [nodes  coords]; 

%%%  Find  connectivity  of  all  triangles  (surface  elements) 
triinfo  =  xmeshinf o (fem, ' out ' , ' elements ' , ' meshtype ' , ' tri ' ) ; 
connectivity  =  triinfo . nodes ' -1 ;  %%%start  with  node  0 

connectivity  =  [connectivity (:, 1)  connectivity (:, 3)  connectivity (:, 6)  ... 

connectivity (:, 2 )  connectivity (:, 5)  connectivity (:,  4 )] ; 

%%%  Eliminate  internal  nodes  from  list,  keeping  surface  nodes. 

AA  =  lel0*ones ( length (allnodecoords ), 4 ) ; 
for  ii  =  1 : length (connectivity) ; 

AA (connectivity (ii, 1) +1, : )  =  allnodecoords (connectivity (ii, 1) +1, :) ; 

AA (connectivity (ii, 2) +1, : )  =  allnodecoords (connectivity (ii, 2) +1, :) ; 

AA (connectivity (ii, 3) +1, : )  =  allnodecoords (connectivity (ii, 3) +1, :); 

AA (connectivity (ii, 4 ) +1,  : )  =  allnodecoords (connectivity (ii, 4 ) +1,  :)  ; 

AA (connectivity (ii, 5) +1, : )  =  allnodecoords (connectivity (ii, 5) +1, : ) ; 

AA (connectivity (ii, 6) +1, : )  =  allnodecoords (connectivity (ii,6)+l,:); 
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end 


AO  =  []  ; 

for  ii  =  1 : length (AA) ; 

if  (AA (ii, 2)  ~=  lelO  ||  AA(ii,3)  ~=  lelO  ||  AA(ii,4)  ~=  lelO); 
AO  =  [AO;  AA (ii,  : ) ]  ; 

end 

end 

nodecoords  =  AO; 

%%%  Save  connectivity  and  coordinate  information  to  file 
fidl  =  f open ( ' kirchhof f_conn . mphtxt ' , ' w ' ) ; 
fid2  =  f open (' kirchhof f_coords . mphtxt w ') ; 

for  ii  =  1 : length (connectivity) ; 

fprintf(fidl,'%d  %d  %d  %d  %d  %d\n ' , [connectivity (ii, : ) ] ) ; 

end 

for  ii  =  1 : length (nodecoords )  ; 

fprintf(fid2,'%d  %f  %f  %f \n ' ,  [node coords (ii,  : ) ] )  ; 

end 

fclose (fidl) ; 
fclose (fid2) ; 


After  the  surface  mesh  files  have  been  generated,  they  can  be  read  to  evaluate  the 
Kirchhoff  integral  using  the  functions  kirchhoff_integral.m,  newell.m,  and  trianglepoint.m, 
written  for  use  in  MATLAB  (ver.  7.2.0.232).  The  inputs  are  the  surface  mesh  file  names  and  the 
wavenumber  of  the  incident  plane  wave. 


oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 


o  o  o 
"6  "o  "o 


kirchhof f_integral .m 
Greg  Dietzen  5/27/08 


ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 


%%%  Calculates  the  backscattering  amplitude  of  an  object  from  an  incident 
%%%  wave  in  the  negative  z  direction  using  the  Kirchhoff  integral.  This 
%%%  integral  is  solved  numerically  from  a  surface  mesh  generated  by 
%%%  kirchhof f_mesh .m.  Inputs  are  the  names  of  the  files  containing 
%%%  connectivity  and  coordinate  information  as  well  as  wavenumber  k. 

function  f  =  kirchhof f_integral (coords_f ilename,  conn_f ilename,  k) ; 

fidl  =  f open (conn_f ilename, ' r ') ; 
fid2  =  f open (coords_f ilename,  ' r ')  ; 

connectivity  =  fscanf (fidl,  ' %d  %d  %d  %d  %d  %d\n ' ,  [6,  inf ] )  '  ; 
coords  =  fscanf (fid2,'%d  %f  %f  %f\n',[4,inf])'; 
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fclose (fidl) ; 
fclose (fid2) ; 

%%%  Determine  the  insonified  region  of  the  surface  mesh  using  Newell's 
%%%  algorithm. 

visible_conn  =  newell (coords,  connectivity); 

kx  =  0 ; 
ky  =  0; 
kz  =  -k; 

k_vector  =  [kx  ky  kz] ; 

lamda  =  2*pi/k; 

%%%  Evaluate  the  integral . 
f  =  0; 

k_unit  =  [kx  ky  kz] / (sqrt (kx^2+ky^2+kz^2) ) ; 
for  ii  =  1 : length (visible_conn) ; 
nodel  =  visible_conn (ii,  1)  ; 
node2  =  visible_conn (ii,  2)  ; 
node3  =  visible_conn (ii, 3)  ; 
xl  =  coords (nodel  +  1, 2)  ; 
yl  =  coords (nodel  +  1 , 3 )  ; 
zl  =  coords (nodel+1, 4 ) ; 
x2  =  coords (node2+l , 2 ) ; 
y2  =  coords (node2+l , 3 ) ; 
z2  =  coords (node2+l , 4 ) ; 
x3  =  coords (node3+l, 2) ; 
y3  =  coords (node3+l , 3 ) ; 
z3  =  coords (node3+l, 4 ) ; 

a  =  sqrt ( (xl-x2)  . A2+ (yl-y2)  . A2+  (zl-z2)  .  A2) ; 
b  =  sqrt ( (x3-x2)  . A2+ (y3-y2)  . A2+  (z3-z2)  .  A2) ; 
c  =  sqrt ( (xl-x3)  . A2+ (yl-y3)  . A2+  (zl-z3)  . A2) ; 
s  =  (a+b+c) . / 2; 

area  =  sqrt (s . * (s-a)  . *  (s-b)  . *  (s-c) ) ; 
vectl2  =  [x2-xl  y2-yl  z2-zl]; 
vectl3  =  [x3-xl  y3-yl  z3-zl]; 

normal  =  [vectl2 ( : , 2 ) . *vectl3 ( : , 3) -vectl2 ( : , 3) . *vectl3 ( : , 2) 
vectl2  ( : , 3)  . *vectl3 ( : , 1 ) -vectl2 ( : , 1)  . *vectl3 ( : , 3)  vectl2 ( : , 1)  . *vectl3 ( : , 2) - 
vectl2 ( : , 2 ) . *vectl3 ( : , 1 ) ] ; 

unit_n  =  normal/sqrt (normal (1) A2+normal (2) A2+normal (3) A2) ; 

n_dot_k  =  -abs (dot (unit_n,  k_unit)); 

rl  =  [xl  yl  zl] ; 

r2  =  [x2  y2  z2]  ; 

r3  =  [x3  y3  z3] ; 

rl_dot_k  =  dot (rl, k_vector) ; 

r2_dot_k  =  dot (r2, k_vector) ; 

r3_dot_k  =  dot (r3, k_vector ) ; 

r_dot_k  =  mean ( [ r l_dot_k  r2_dot_k  r3_dot_k] ) ; 
fO  =  i/lamda*n_dot_k*exp (i*2*r_dot_k) *area; 
f  =  f+fO; 

end 
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%%%  Determines  the  insonified  region  of  a  surface  mesh  using  a  modified 
%%%  Newell's  algorithm.  Incoming  wave  is  in  the  negative  z  direction. 

function  visible_conn  =  newell (coordinates,  connectivity); 

warning  off  MATLAB : divideByZero 

%%%  Number  triangular  elements  for  labeling  purposes, 
connectivity  =  [[ 0 : length (connectivity) -1 ] '  connectivity]; 

%%%  Sample  size  for  rasterization  test, 
samplesize  =  20; 

%%%  Determine  min  and  max  x,y,  and  z  values  for  each  element  and  sort  by 
%%%  max  z . 

maxmin  =  zeros  (length (connectivity) , 7) ; 
for  ii  =  1 : length (connectivity) ; 
nodel  =  connectivity (ii,  2)  ; 
node2  =  connectivity (ii, 3) ; 
node3  =  connectivity (ii, 4 )  ; 

maxx  =  max ( [coordinates (nodel+1, 2)  coordinates (node2+l, 2)  ... 

coordinates (node3  +  l,  2) ]  )  ; 

maxy  =  max ( [coordinates (nodel  +  1,  3)  coordinates (node2  +  l , 3 )  ... 

coordinates (node3  +  l,  3)  ]  )  ; 

maxz  =  max ( [coordinates (nodel+1, 4 )  coordinates (node2+l , 4 )  ... 

coordinates (node 3  +  1 ,  4 )  ]  )  ; 

minx  =  min ( [coordinates (nodel+1, 2)  coordinates (node2+l , 2 )  ... 

coordinates (node3+l, 2) ] ) ; 

miny  =  min ([ coordinates (nodel+1 , 3 )  coordinates (node2+l, 3)  ... 

coordinates (node3+l, 3) ] ) ; 

minz  =  min ( [coordinates (nodel+1, 4 )  coordinates (node2+l, 4 )  ... 

coordinates (node3+l, 4) ] ) ; 

maxmin (ii,:)  =  [ii-1  minx  maxx  miny  maxy  minz  maxz]; 

end 

sortmaxmin  =  sortrows (maxmin, 6 ) ; 
sortmaxmin  =  flipud ( sortrows ( sortmaxmin, 7 )) ; 

%%%  Minmax  tests  in  x  and  y  and  rasterization  test, 
ii  =  1; 

while  ii  <=  length (sortmaxmin) ; 
maxx  =  sortmaxmin (ii,  3)  ; 
minx  =  sortmaxmin (ii, 2) ; 
maxy  =  sortmaxmin (ii, 5) ; 
miny  =  sortmaxmin (ii, 4 ) ; 
jj  =  ii+1 ; 

while  jj  <=  length (sortmaxmin) ; 
maxx2  =  sortmaxmin ( j j , 3 )  ; 
minx2  =  sortmaxmin ( j j  ,  2 )  ; 
if  (minx  >=  maxx2  | |  maxx  <=  minx2) ; 

jj  =  jj+i; 

else 

maxy2  =  sortmaxmin ( j j , 5 ) ; 
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miny2  =  sortmaxmin ( j j , 4 ) ; 

if  (miny  >=  maxy2  | |  maxy  <=  miny2) ; 

jj  =  jj+i; 

else 

nodel  =  connectivity (sortmaxmin (ii, 1) +1, 2) ; 
node2  =  connectivity (sortmaxmin (ii, 1) +1, 3) ; 
node3  =  connect ivity ( sortmaxmin ( ii, 1 ) +1 , 4 ) ; 
node4  =  connect ivity ( sortmaxmin ( j j , 1 ) +1 , 2 ) ; 
node5  =  connect ivity ( sortmaxmin ( j j , 1 ) +1 , 3 ) ; 
node6  =  connect ivity ( sortmaxmin ( j j , 1 ) +1 , 4 ) ; 
raster  =  zeros (samplesize, samplesize)  ; 
xxx  =  [minx+ (maxx-minx) /samplesize : (maxx-minx) /.. . 
samplesize :maxx] ; 

yyy  =  fliplr ( [miny+ (maxx-minx) /samplesize : (maxy-miny) / . 

samplesize :maxy] ) ; 
for  mm  =  1 : length  (xxx) ; 

for  nn  =  1 : length (yyy ) ; 
x  =  xxx (mm) ; 

y  =  yyy (nn) ; 

xl  =  coordinates (nodel+1, 2) ; 
yl  =  coordinates (nodel+1, 3) ; 
x2  =  coordinates (node2+l , 2 ) ; 
y2  =  coordinates (node2+l , 3 ) ; 
x3  =  coordinates (node3+l, 2) ; 
y3  =  coordinates (node3+l , 3 ) ; 

tril  =  trianglepoint (xl, yl, x2, y2, x3, y3, x, y) ; 
raster (nn, mm)  =  tril; 

end 

end 

overlap  =  0; 

for  mm  =  1 : length (xxx) ; 

for  nn  =  1 : length (yyy ) ; 
x  =  xxx (mm) ; 
y  =  yyy (nn) ; 

x4  =  coordinates (node4+l, 2) ; 
y4  =  coordinates (node4+l, 3) ; 
x5  =  coordinates (node5+l , 2 ) ; 
y5  =  coordinates (node5+l , 3 ) ; 
x 6  =  coordinates (node6+l, 2) ; 
y6  =  coordinates (node6+l, 3) ; 

tri2  =  trianglepoint (x4 , y4 , x5, y5, x6, y6, x, y) ; 
if  (tri2  ==  1  &&  raster (nn, mm)  ==  1); 
overlap  =  1; 

end 

end 

end 

if  overlap  ==  1; 

sortmaxmin ( j j ,: )  =  []; 

else 

jj  =  jj+i; 

end 

end 

end 

end 

ii  =  ii+1 ; 

end 
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visible_conn  =  zeros (length (sortmaxmin) , 6) ; 

for  ii  =  1 : length (sortmaxmin)  ; 

visible_conn (ii, : )  =  connectivity (sortmaxmin (ii, 1) +1, 2 : 7) ; 

end 


%%%  trianglepoint . m 
%%%  Greg  Dietzen  5/27/08 


oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 


%%%  This  function  determines  whether  a  point  lies  within  a  triangle.  Both 
%%%  the  point  (defined  by  inputs  x  and  y)  and  the  triangle  (defined  by 
%%%  vertices  (xl,yl), (x2,y2),  and  (x3,y3))  are  in  the  x  y  plane.  The  output 
%%%  is  1  if  the  point  is  within  the  triangle  and  0  if  it  is  not.  If  the 
%%%  point  lies  on  the  edge  of  the  triangle,  it  is  not  considered  part  of 


function  ans  =  trianglepoint (xl, yl, x2, y2, x3, y3, x, y) ; 
ml2  =  (y2-yl) / (x2-xl) ; 

if  (isnan(ml2)  ==  1  | |  ml2  ==  Inf  | |  ml2  ==  -Inf) ; 
if  x3  >  xl; 

if  x  >  xl ; 

m23  =  (y3-y2 ) / (x3-x2 ) ; 

if  (isnan(m23)  ==  1  | |  m23  ==  Inf  | |  m23  ==  -Inf) ; 
if  xl  >  x2; 

if  x  >  x2; 

m31  =  (yl-y3) / (xl-x3) ; 

if  (isnan(m31)  ==  1  | |  m31  ==  Inf  | |  m31  ==  -Inf) ; 
if  x2  >  x3; 

if  x  >  x3; 
ans  =  1; 

else 

ans  =  0; 

end 

else 

if  x  <  x3; 
ans  =  1; 

else 

ans  =  0; 

end 

end 

else 

yy2  =  m31* (x2-x3) +y3; 
if  yy2  <  y2; 

yy  =  m31* (x-x3) +y3; 
if  yy  <  y; 

ans  =  1; 

else 

ans  =  0; 

end 

else 

yy  =  m31 * (x-x3 ) +y3; 
if  yy  >  y;  %%%point  is  below  3-1 
ans  =  1; 
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else 


0; 


ans  = 

end 

end 

end 

else 

ans  =  0; 

end 

else 

if  x  <  x2; 

m31  =  (yl-y3) / (xl-x3) ; 
if  (isnan(m31)  ==  1  | |  m31  ==  Inf  | |  m31  ==  -Inf) 
if  x2  >  x3; 

if  x  >  x3; 
ans  =  1; 

else 

ans  =  0; 

end 

else 

if  x  <  x3; 
ans  =  1; 

else 

ans  =  0; 

end 

end 

else 

yy2  =  m31* (x2-x3) +y3; 
if  yy2  <  y2; 

yy  =  m31* (x-x3) +y3; 
if  yy  <  y; 

ans  =  1; 

else 

ans  =  0; 

end 

else 

yy  =  m31* (x-x3) +y3; 
if  yy  >  y; 

ans  =  1; 

else 

ans  =  0; 

end 

end 

end 

else 

ans  =  0; 

end 

end 

else 

yyl  =  m23* (xl-x2 ) +y2 ; 
if  yyl  <  yl; 

yy  =  m23* (x-x2 ) +y2 ; 
if  yy  <  y; 

m31  =  (yl-y3 ) / (xl-x3 )  ; 

if  (isnan(m31)  ==  1  |  f  m31  ==  Inf  | |  m31  ==  -Inf) 
if  x2  >  x3; 

if  x  >  x3; 
ans  =  1; 
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else 


0; 


ans  = 

end 

else 

if  x  <  x3; 
ans  =  1; 

else 

ans  =  0; 

end 

end 

else 

yy2  =  m31* (x2-x3) +y3; 
if  yy2  <  y2; 

yy  =  m31* (x-x3) +y3; 

if  yy  <  y; 

ans  =  1; 

else 

ans  =  0; 

end 

else 

yy  =  m31* (x-x3) +y3; 
if  yy  >  y; 

ans  =  1; 

else 

ans  =  0; 

end 

end 

end 

else 

ans  =  0; 

end 

else 

yy  =  m23* (x-x2 ) +y2; 
if  yy  >  y; 

m31  =  (yl-y3 ) / (xl-x3 )  ; 

if  (isnan(m31)  ==  1  | |  m31  ==  Inf  | |  m31  ==  -Inf) 
if  x2  >  x3; 

if  x  >  x3; 
ans  =  1; 

else 

ans  =  0; 

end 

else 

if  x  <  x3; 
ans  =  1; 

else 

ans  =  0; 

end 

end 

else 

yy2  =  m31* (x2-x3) +y3; 
if  yy2  <  y2; 

yy  =  m31* (x-x3) +y3; 

if  yy  <  y; 

ans  =  1; 

else 

ans  =  0; 
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end 


else 

yy  =  m31* (x-x3) +y3; 
if  yy  >  y; 

ans  =  1; 

else 

ans  =  0; 

end 

end 

end 

else 

ans  =  0; 

end 

end 

end 

else 

ans  =  0; 

end 

else 

if  x  <  xl ; 

m23  =  (y3-y2 ) / (x3-x2 ) ; 

if  (isnan(m23)  ==  1  | |  m23  ==  Inf  | |  m23  ==  -Inf) ; 
if  xl  >  x2; 

if  x  >  x2; 

m31  =  (yl-y3 ) / (xl-x3 ) ; 

if  (isnan(m31)  ==  1  | |  m31  ==  Inf  | |  m31  ==  -Inf) 
if  x2  >  x3; 

if  x  >  x3; 
ans  =  1; 

else 

ans  =  0; 

end 

else 

if  x  <  x3; 
ans  =  1; 

else 

ans  =  0; 

end 

end 

else 

yy2  =  m31* (x2-x3) +y3; 
if  yy2  <  y2; 

yy  =  m31* (x-x3) +y3; 
if  yy  <  y; 

ans  =  1; 

else 

ans  =  0; 

end 

else 

yy  =  m31 * (x-x3 ) +y3; 
if  yy  >  y; 

ans  =  1; 

else 

ans  =  0; 

end 

end 

end 
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else 


0; 


ans  = 

end 

else 

if  x  <  x2; 

m31  =  (yl-y3) / (xl-x3) ; 

if  (isnan(m31)  ==  1  | |  m31  ==  Inf  | |  m31  ==  -Inf) 
if  x2  >  x3; 

if  x  >  x3; 
ans  =  1; 

else 

ans  =  0; 

end 

else 

if  x  <  x3; 
ans  =  1; 

else 

ans  =  0; 

end 

end 

else 

yy2  =  m31* (x2-x3) +y3; 
if  yy2  <  y2; 

yy  =  m31* (x-x3) +y3; 
if  yy  <  y; 

ans  =  1; 

else 

ans  =  0; 

end 

else 

yy  =  m31* (x-x3) +y3; 
if  yy  >  y; 

ans  =  1; 

else 

ans  =  0; 

end 

end 

end 

else 

ans  =  0; 

end 

end 

else 

yyl  =  m23* (xl-x2 ) +y2 ; 
if  yyl  <  yl; 

yy  =  m23* (x-x2 ) +y2 ; 
if  yy  <  y; 

m31  =  (yl-y3 ) / (xl-x3 )  ; 

if  (isnan(m31)  ==  1  ]}  m31  ==  Inf  | |  m31  ==  -Inf) 
if  x2  >  x3; 

if  x  >  x3; 
ans  =  1; 

else 

ans  =  0; 

end 

else 

if  x  <  x3; 
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ans 


1; 


else 

ans  =  0; 

end 

end 

else 

yy2  =  m31* (x2-x3) +y3; 
if  yy2  <  y2; 

yy  =  m31* (x-x3) +y3; 

if  yy  <  y; 

ans  =  1; 

else 

ans  =  0; 

end 

else 

yy  =  m31* (x-x3) +y3; 
if  yy  >  y; 

ans  =  1; 

else 

ans  =  0; 

end 

end 

end 

else 

ans  =  0; 

end 

else 

yy  =  m23* (x-x2 ) +y2; 
if  yy  >  y; 

m31  =  (yl-y3 ) / (xl-x3 )  ; 

if  (isnan(m31)  ==  1  | |  m31  ==  Inf  | |  m31  ==  -Inf) 
if  x2  >  x3; 

if  x  >  x3; 
ans  =  1; 

else 

ans  =  0; 

end 

else 

if  x  <  x3; 
ans  =  1; 

else 

ans  =  0; 

end 

end 

else 

yy2  =  m31* (x2-x3) +y3; 
if  yy2  <  y2; 

yy  =  m31* (x-x3) +y3; 

if  yy  <  y; 

ans  =  1; 

else 

ans  =  0; 

end 

else 

yy  =  m31* (x-x3) +y3; 
if  yy  >  y; 

ans  =  1; 
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else 


0; 


ans  = 


end 

end 

end 

else 

ans  =  0; 

end 

end 

end 

else 

ans  =  0; 

end 

end 

else 

yy3  =  ml2* (x3-xl ) +yl ; 
if  yy3  <  y3; 

yy  =  ml2* (x-xl) +yl; 
if  yy  <  y; 

m23  =  (y3-y2 ) / (x3-x2 )  ; 

if  (isnan(m23)  ==  1  | |  m23  ==  Inf  | |  m23 
if  xl  >  x2; 

if  x  >  x2; 

m31  =  (yl-y3) / (xl-x3) ; 
if  (isnan(m31)  ==  1  | |  m31  = 
if  x2  >  x3; 

if  x  >  x3; 
ans  =  1; 

else 

ans  =  0; 

end 

else 

if  x  <  x3; 
ans  =  1; 

else 

ans  =  0; 

end 


=  -Inf)  ; 


Inf  |  |  m31  =: 


end 

else 

yy2  =  m31* (x2-x3) +y3; 
if  yy2  <  y2; 

yy  =  m31* (x-x3) +y3; 
if  yy  <  y; 

ans  =  1; 

else 

ans  =  0; 

end 

else 

yy  =  m31* (x-x3) +y3; 
if  yy  >  y; 

ans  =  1; 

else 

ans  =  0; 

end 

end 

end 

else 


-Inf) 
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end 


ans 


0; 


else 

if  x  <  x2; 

m31  =  (yl-y3) / (xl-x3) ; 

if  (isnan(m31)  ==  1  | |  m31  ==  Inf  | |  m31  ==  -Inf) 
if  x2  >  x3; 

if  x  >  x3; 
ans  =  1; 

else 

ans  =  0; 

end 

else 

if  x  <  x3; 
ans  =  1; 

else 

ans  =  0; 

end 

end 

else 

yy2  =  m31* (x2-x3) +y3; 
if  yy2  <  y2; 

yy  =  m31* (x-x3) +y3; 
if  yy  <  y; 

ans  =  1; 

else 

ans  =  0; 

end 

else 

yy  =  m31* (x-x3) +y3; 
if  yy  >  y; 

ans  =  1; 

else 

ans  =  0; 

end 

end 

end 

else 

ans  =  0; 

end 

end 

else 

yyl  =  m23* (xl-x2 ) +y2 ; 
if  yyl  <  yl; 

yy  =  m23* (x-x2 ) +y2; 
if  yy  <  y; 

m31  =  (yl-y3 ) / (xl-x3 )  ; 

if  (isnan(m31)  ==  1  | |  m31  ==  Inf  | |  m31  ==  -Inf) 
if  x2  >  x3; 

if  x  >  x3; 
ans  =  1; 

else 

ans  =  0; 

end 

else 

if  x  <  x3; 
ans  =  1; 


149 


else 


0; 


ans  = 

end 

end 

else 

yy2  =  m31* (x2-x3) +y3; 
if  yy2  <  y2; 

yy  =  m31* (x-x3) +y3; 

if  yy  <  y; 

ans  =  1; 

else 

ans  =  0; 

end 

else 

yy  =  m31* (x-x3) +y3; 
if  yy  >  y; 

ans  =  1; 

else 

ans  =  0; 

end 

end 

end 

else 

ans  =  0; 

end 

else 

yy  =  m23* (x-x2 ) +y2; 
if  yy  >  y; 

m31  =  (yl-y3 ) / (xl-x3 )  ; 

if  (isnan(m31)  ==  1  | |  m31  ==  Inf  | |  m31  ==  -Inf) 
if  x2  >  x3; 

if  x  >  x3; 
ans  =  1; 

else 

ans  =  0; 

end 

else 

if  x  <  x3; 
ans  =  1; 

else 

ans  =  0; 

end 

end 

else 

yy2  =  m31* (x2-x3) +y3; 
if  yy2  <  y2; 

yy  =  m31* (x-x3) +y3; 
if  yy  <  y; 

ans  =  1; 

else 

ans  =  0; 

end 

else 

yy  =  m31* (x-x3) +y3; 
if  yy  >  y; 

ans  =  1; 

else 
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end 


ans 


0; 


end 

end 

else 

ans  =  0; 

end 

end 

end 

else 

ans  =  0; 

end 

else 

yy  =  ml2* (x-xl) +yl; 
if  yy  >  y; 

m23  =  (y3-y2 ) / (x3-x2 )  ; 

if  (isnan(m23)  ==  1  |  J  m23  ==  Inf  |  |  m23  ==  -Inf)  ; 
if  xl  >  x2; 

if  x  >  x2; 

m31  =  (yl-y3) / (xl-x3) ; 

if  (isnan(m31)  ==  1  | |  m31  ==  Inf  | |  m31  ==  -Inf) 
if  x2  >  x3; 

if  x  >  x3; 
ans  =  1; 

else 

ans  =  0; 

end 

else 

if  x  <  x3 ; 
ans  =  1; 

else 

ans  =  0; 

end 

end 

else 

yy2  =  m31* (x2-x3) +y3; 
if  yy2  <  y2; 

yy  =  m31* (x-x3) +y3; 
if  yy  <  y; 

ans  =  1; 

else 

ans  =  0; 

end 

else 

yy  =  m31 * (x-x3 ) +y3; 
if  yy  >  y; 

ans  =  1; 

else 

ans  =  0; 

end 

end 

end 

else 

ans  =  0; 

end 

else 

if  x  <  x2; 
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m31  =  (yl-y3) / (xl-x3) ; 

if  (isnan(m31)  ==  1  |  j  m31  ==  Inf  | |  m31  ==  -Inf) 
if  x2  >  x3; 

if  x  >  x3; 
ans  =  1; 

else 

ans  =  0; 

end 

else 

if  x  <  x3; 
ans  =  1; 

else 

ans  =  0; 

end 

end 

else 

yy2  =  m31* (x2-x3) +y3; 
if  yy2  <  y2; 

yy  =  m31* (x-x3) +y3; 
if  yy  <  y; 

ans  =  1; 

else 

ans  =  0; 

end 

else 

yy  =  m31* (x-x3) +y3; 
if  yy  >  y; 

ans  =  1; 

else 

ans  =  0; 

end 

end 

end 

else 

ans  =  0; 

end 

end 

else 

yyl  =  m23* (xl-x2 ) +y2 ; 
if  yyl  <  yl; 

yy  =  m23* (x-x2 ) +y2 ; 
if  yy  <  y; 

m31  =  (yl-y3 ) / (xl-x3 )  ; 

if  (isnan(m31)  ==  1  | |  m31  ==  Inf  | |  m31  ==  -Inf) 
if  x2  >  x3; 

if  x  >  x3; 
ans  =  1; 

else 

ans  =  0; 

end 

else 

if  x  <  x3; 
ans  =  1; 

else 

ans  =  0; 

end 

end 
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else 


yy2  =  m31* (x2-x3) +y3; 
if  yy2  <  y2; 

yy  =  m31* (x-x3) +y3; 

if  yy  <  y; 

ans  =  1; 

else 

ans  =  0; 

end 

else 

yy  =  m31* (x-x3) +y3; 
if  yy  >  y; 

ans  =  1; 

else 

ans  =  0; 

end 

end 

end 

else 

ans  =  0; 

end 

else 

yy  =  m23* (x-x2 ) +y2; 
if  yy  >  y; 

m31  =  (yl-y3 ) / (xl-x3 )  ; 

if  (isnan(m31)  ==  1  | |  m31  ==  Inf  | |  m31  ==  -Inf) 
if  x2  >  x3; 

if  x  >  x3; 
ans  =  1; 

else 

ans  =  0; 

end 

else 

if  x  <  x3; 
ans  =  1; 

else 

ans  =  0; 

end 

end 

else 

yy2  =  m31* (x2-x3) +y3; 
if  yy2  <  y2; 

yy  =  m31* (x-x3) +y3; 
if  yy  <  y; 

ans  =  1; 

else 

ans  =  0; 

end 

else 

yy  =  m31* (x-x3) +y3; 
if  yy  >  y; 

ans  =  1; 

else 

ans  =  0; 

end 

end 

end 
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else 


ans  =  0; 

end 

end 

end 

else 

ans  =  0; 

end 

end 

end 
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